! modal_aero_coag.F90


!----------------------------------------------------------------------
!BOP
!
! !MODULE: modal_aero_coag --- modal aerosol coagulation
!
! !INTERFACE:
   module modal_aero_coag

! !USES:
   use shr_kind_mod,    only:  r8 => shr_kind_r8
   use shr_kind_mod,    only:  r4 => shr_kind_r4
   use chem_mods,       only:  gas_pcnst
   use modal_aero_data, only:  maxd_aspectype

  implicit none
  private
  save

! !PUBLIC MEMBER FUNCTIONS:
  public modal_aero_coag_sub, modal_aero_coag_init

! !PUBLIC DATA MEMBERS:
  integer, parameter :: pcnstxx = gas_pcnst

#if ( defined MODAL_AERO_7MODE )
  integer, parameter, public :: pair_option_acoag = 3
#elif ( defined MODAL_AERO_3MODE )
  integer, parameter, public :: pair_option_acoag = 1
#endif
! specifies pairs of modes for which coagulation is calculated
!   1 -- [aitken-->accum] 
!   2 -- [aitken-->accum], and [pcarbon-->accum]
!   3 -- [aitken-->accum], [pcarbon-->accum], 
!        and [aitken-->pcarbon--(aging)-->accum] 
! other -- do no coag

  integer, parameter, public :: maxpair_acoag = 10
  integer, parameter, public :: maxspec_acoag = maxd_aspectype

  integer, public :: npair_acoag
  integer, public :: modefrm_acoag(maxpair_acoag)
  integer, public :: modetoo_acoag(maxpair_acoag)
  integer, public :: modetooeff_acoag(maxpair_acoag)
  integer, public :: nspecfrm_acoag(maxpair_acoag)
  integer, public :: lspecfrm_acoag(maxspec_acoag,maxpair_acoag)
  integer, public :: lspectoo_acoag(maxspec_acoag,maxpair_acoag)

! !DESCRIPTION: This module implements ...
!
! !REVISION HISTORY:
!
!   RCE 07.04.13:  Adapted from MIRAGE2 code
!
!EOP
!----------------------------------------------------------------------
!BOC

! list private module data here

!EOC
!----------------------------------------------------------------------


  contains
                                                                                                                                            
!----------------------------------------------------------------------
!BOP
! !ROUTINE:  modal_aero_coag_sub --- ...
!
! !INTERFACE:
   subroutine modal_aero_coag_sub(                               &
                        lchnk,    ncol,     nstep,               &
                        loffset,  deltat_main,                   &
                        latndx,   lonndx,                        &
                        t,        pmid,     pdel,                &
                        q,                                       &
                        dgncur_a,           dgncur_awet,         &
                        wetdens_a                                )


!----------------------------------------------------------------------
!   Authors: R. Easter
!----------------------------------------------------------------------

! !USES:
   use mo_constants,     only: pi
   use modal_aero_data
   use modal_aero_gasaerexch, only:  n_so4_monolayers_pcage, &
                                     soa_equivso4_factor

   use abortutils,       only: endrun
   use cam_history,      only: outfld, fieldname_len
   use chem_mods,        only: adv_mass
   use constituents,     only: pcnst, cnst_name
   use physconst,        only: gravit, mwdry, r_universal
   use ppgrid,           only: pcols, pver
   use spmd_utils,       only: iam, masterproc

   implicit none

! !PARAMETERS:
   integer, intent(in)  :: lchnk            ! chunk identifier
   integer, intent(in)  :: ncol             ! number of columns in chunk
   integer, intent(in)  :: nstep            ! model step
   integer, intent(in)  :: loffset          ! offset applied to modal aero "pointers"
   integer, intent(in)  :: latndx(pcols), lonndx(pcols) 

   real(r8), intent(in) :: deltat_main      ! model timestep (s)

   real(r8), intent(in) :: t(pcols,pver)    ! temperature (K)
   real(r8), intent(in) :: pmid(pcols,pver) ! pressure at model levels (Pa)
   real(r8), intent(in) :: pdel(pcols,pver) ! pressure thickness of levels (Pa)

   real(r8), intent(inout) :: q(ncol,pver,pcnstxx) 
                                            ! tracer mixing ratio (TMR) array
                                            ! *** MUST BE mol/mol-air or #/mol-air
                                            ! *** NOTE ncol & pcnstxx dimensions
   real(r8), intent(in) :: dgncur_a(pcols,pver,ntot_amode)
                                 ! dry geo. mean dia. (m) of number distrib.
   real(r8), intent(in) :: dgncur_awet(pcols,pver,ntot_amode)
                                 ! wet geo. mean dia. (m) of number distrib.
   real(r8), intent(in) :: wetdens_a(pcols,pver,ntot_amode) 
                                 ! density of wet aerosol (kg/m3)

! !DESCRIPTION: 
!   computes changes due to coagulation involving
!	aitken mode (modeptr_aitken) with accumulation mode (modeptr_accum)
!   this version will 
!	compute changes to mass and number, but not to surface area
!	calculates coagulation rate coefficients using either
!	    new CMAQ V4.6 fast method
!	    older cmaq slow method (direct gauss-hermite quadrature)
!
! !REVISION HISTORY:
!   RCE 07.04.15:  Adapted from MIRAGE2 code and CMAQ V4.6 code
!
!EOP
!----------------------------------------------------------------------
!BOC

! local variables
	integer :: i, iok, ipair, ip_aitacc, ip_aitpca, ip_pcaacc, iq
	integer :: idomode(ntot_amode), iselfcoagdone(ntot_amode)
	integer :: jfreqcoag
	integer :: k
	integer :: l, l1, l2, la, lmz, lsfrm, lstoo, lunout
	integer :: modefrm, modetoo, mait, macc, mpca
	integer ::  n, nfreqcoag


	integer, save :: nerr = 0       ! number of errors for entire run
	integer, save :: nerrmax = 9999 ! maximum number of errors before abort
	integer, parameter :: ldiag1=-1, ldiag2=-1, ldiag3=-1

	logical, parameter :: fastcoag_flag = .true. ! selects coag rate-coef method

	real(r8) :: aircon
      	real(r8) :: deltat, deltatinv_main
      	real(r8) :: dr_so4_monolayers_pcage
	real(r8) :: dryvol_a(pcols,pver,ntot_amode)
      	real(r8) :: dumexp, dumloss, dumprod
	real(r8) :: dumsfc_frm_old, dumsfc_frm_new
	real(r8) :: dum_m2v
	real(r8) :: fac_m2v_aitage(maxd_aspectype), fac_m2v_pcarbon(maxd_aspectype)
	real(r8) :: fac_volsfc_pcarbon
	real(r8) :: lnsg_frm, lnsg_too
	real(r8) :: sg_frm, sg_too
	real(r8) :: tmpa, tmpb, tmpc, tmpf, tmpg, tmph, tmpn
	real(r8) :: tmp1, tmp2
	real(r8) :: tmp_qold
	real(r8) :: v2ncur_a_tmp
	real(r8) :: vol_core, vol_shell
	real(r8) :: wetdens_frm, wetdens_too, wetdgnum_frm, wetdgnum_too
	real(r8) :: xbetaij0, xbetaij2i, xbetaij2j, xbetaij3, &
                    xbetaii0, xbetaii2,  xbetajj0,  xbetajj2     
      	real(r8) :: xferamt, xferfracvol, xferfrac_pcage, xferfrac_max
	real(r8) :: xnumbconc(ntot_amode)
	real(r8) :: xnumbconcavg(ntot_amode), xnumbconcnew(ntot_amode)
	real(r8) :: ybetaij0(maxpair_acoag), ybetaij3(maxpair_acoag)
	real(r8) :: ybetaii0(maxpair_acoag), ybetajj0(maxpair_acoag)

        real(r8) :: dqdt(ncol,pver,pcnstxx)  ! TMR "dq/dt" array - NOTE dims
        logical  :: dotend(pcnst)            ! identifies the species that
                                             ! tendencies are computed for
	real(r8) :: qsrflx(pcols)

        character(len=fieldname_len)   :: tmpname
        character(len=fieldname_len+3) :: fieldname

! begin
!   check if any coagulation pairs exist
	if (npair_acoag <= 0) return

!--------------------------------------------------------------------------------
   if (ldiag1 > 0) then
   if (nstep <= 3) then
   do i = 1, ncol
   if (lonndx(i) /= 37) cycle
   if (latndx(i) /= 23) cycle
   if (nstep > 3)       cycle
   write( *, '(/a,i7,i5,2(2x,2i5))' )   &
         '*** modal_aero_coag_sub -- nstep, iam, lat, lon, pcols, ncol =',   &
         nstep, iam, latndx(i), lonndx(i), pcols, ncol
   end do
   end if
!  if (ncol /= -999888777) return
   if (nstep > 3) call endrun( 'modal_aero_coag_sub -- nstep>3 testing halt' )
   end if   ! (ldiag1 > 0)
!--------------------------------------------------------------------------------

	dotend(:) = .false.
	dqdt(1:ncol,:,:) = 0.0

	lunout = 6


!
!   determine if coagulation will be done on this time-step
!   currently coagulation is done every 3 hours
!
!	deltat = 3600.0*3.0
	deltat = deltat_main
	nfreqcoag = max( 1, nint( deltat/deltat_main ) )
	jfreqcoag = nfreqcoag/2
	xferfrac_max = 1.0_r8 - 10.0_r8*epsilon(1.0_r8)   ! 1-eps

	if (nfreqcoag .gt. 1) then
	    if ( mod(nstep,nfreqcoag) .ne. jfreqcoag ) return
	end if

!
!   set idomode
!
	idomode(:) = 0
	do ipair = 1, npair_acoag
	    idomode(modefrm_acoag(ipair)) = 1
	    idomode(modetoo_acoag(ipair)) = 1
	end do

!
!   other init
!
	macc = modeptr_accum
	mait = modeptr_aitken
	mpca = modeptr_pcarbon

	fac_m2v_aitage(:) = 0.0
	fac_m2v_pcarbon(:) = 0.0
	if (pair_option_acoag == 3) then
!   following ipair definitions MUST BE CONSISTENT with
!   the coding in modal_aero_coag_init for pair_option_acoag == 3
	    ip_aitacc = 1
	    ip_pcaacc = 2
	    ip_aitpca = 3

	! use 1 mol (bi-)sulfate = 65 cm^3 --> 1 molecule = (4.76e-10 m)^3
	    dr_so4_monolayers_pcage = n_so4_monolayers_pcage * 4.76e-10

	    ipair = ip_aitpca
	    do iq = 1, nspecfrm_acoag(ipair)
		lsfrm = lspecfrm_acoag(iq,ipair)
		if (lsfrm == lptr_so4_a_amode(mait)) then
		    fac_m2v_aitage(iq) = specmw_so4_amode / specdens_so4_amode
		else if (lsfrm == lptr_nh4_a_amode(mait)) then
		    fac_m2v_aitage(iq) = specmw_nh4_amode / specdens_nh4_amode
		else if (lsfrm == lptr_soa_a_amode(mait)) then
		    fac_m2v_aitage(iq) = soa_equivso4_factor*   &
                                        (specmw_soa_amode / specdens_soa_amode)
!   for soa, the soa_equivso4_factor converts the soa volume into an
!	so4(+nh4) volume that has same hygroscopicity contribution as soa
!   this allows aging calculations to be done in terms of the amount
!	of (equivalent) so4(+nh4) in the shell
!   (see modal_aero_gasaerexch)
		end if
	    end do
	    
	    do l = 1, nspec_amode(mpca)
		l2 = lspectype_amode(l,mpca)
!   fac_m2v converts (kmol-AP/kmol-air) to (m3-AP/kmol-air)
!		[m3-AP/kmol-AP]    = [kg-AP/kmol-AP]  / [kg-AP/m3-AP]
		fac_m2v_pcarbon(l) = specmw_amode(l2) / specdens_amode(l2)
	    end do

	    fac_volsfc_pcarbon = exp( 2.5*(alnsg_amode(mpca)**2) )
	else
	    ip_aitacc = -999888777
	    ip_pcaacc = -999888777
	    ip_aitpca = -999888777
	end if

!
!   loop over levels and columns to calc the coagulation
!
!   integrate coagulation changes over deltat = nfreqcoag*deltat_main
!   then compute tendencies as
!      dqdt = (q(t+deltat) - q(t))/deltat_main
!   because tendencies are applied (in physics_update) over deltat_main
!
	deltat = nfreqcoag*deltat_main
	deltatinv_main = 1.0_r8/(deltat_main*(1.0_r8 + 1.0e-15_r8))

main_k: do k = 1, pver
main_i: do i = 1, ncol

!   air molar density (kmol/m3)
	aircon = (pmid(i,k)/(r_universal*t(i,k)))

!   calculate number conc. (#/m3) for modes doing coagulation
	do n = 1, ntot_amode
	    if (idomode(n) .gt. 0) then
		xnumbconc(n) = q(i,k,numptr_amode(n)-loffset)*aircon
		xnumbconc(n) = max( 0.0_r8, xnumbconc(n) ) 
	    end if
	    iselfcoagdone(n) = 0
	end do

!
!   calculate coagulation rates for each pair
!
main_ipair1: do ipair = 1, npair_acoag

	modefrm = modefrm_acoag(ipair)
	modetoo = modetoo_acoag(ipair)

!
! compute coagulation rates using cmaq "fast" method
!    (based on E. Whitby's approximation approach)
! here subr. arguments are all in mks unit
!
        call getcoags_wrapper_f(                                       &
          t(i,k), pmid(i,k),                                           &
          dgncur_awet(i,k,modefrm),     dgncur_awet(i,k,modetoo),      &
          sigmag_amode(modefrm),        sigmag_amode(modetoo),         &
          alnsg_amode(modefrm),         alnsg_amode(modetoo),          &
          wetdens_a(i,k,modefrm),       wetdens_a(i,k,modetoo),        &
          xbetaij0, xbetaij2i, xbetaij2j, xbetaij3,                    &
          xbetaii0, xbetaii2,  xbetajj0,  xbetajj2                     )


!   test diagnostics begin --------------------------------------------
 	if (ldiag2 > 0) then
 	if (nstep <= 3) then
 	if ((lonndx(i) == 37) .and. (latndx(i) == 23)) then
 	if ((mod(k-1,5) == 0) .or. (k>=23)) then

	wetdgnum_frm = dgncur_awet(i,k,modefrm)
	wetdgnum_too = dgncur_awet(i,k,modetoo)
	wetdens_frm  = wetdens_a(i,k,modefrm)
	wetdens_too  = wetdens_a(i,k,modetoo)
	sg_frm   = sigmag_amode(modefrm)
	sg_too   = sigmag_amode(modetoo)
	lnsg_frm = alnsg_amode(modefrm)
	lnsg_too = alnsg_amode(modetoo)

        call getcoags_wrapper_f(                                       &
          t(i,k), pmid(i,k),                                           &
          wetdgnum_frm,                   wetdgnum_too,                &
          sg_frm,                         sg_too,                      &
          lnsg_frm,                       lnsg_too,                    &
          wetdens_frm,                    wetdens_too,                 &
          xbetaij0, xbetaij2i, xbetaij2j, xbetaij3,                    &
          xbetaii0, xbetaii2,  xbetajj0,  xbetajj2                     )


 	    write(lunout,9801)
 	    write(lunout,9810) 'nstep,lat,lon,k,ipair   ',   &
 		nstep, latndx(i), lonndx(i), k, ipair
 	    write(lunout,9820) 'tk, pmb, aircon, pdel   ',   &
 		t(i,k), pmid(i,k)*1.0e-2, aircon, pdel(i,k)*1.0e-2
 	    write(lunout,9820) 'wetdens-cgs, sg      f/t',   &
 		wetdens_frm*1.0e-3, wetdens_too*1.0e-3,   &
 		sg_frm, sg_too
 	    write(lunout,9820) 'dgnwet-um, dgndry-um f/t',   &
 		1.0e6*wetdgnum_frm, 1.0e6*wetdgnum_too,   &
 		1.0e6*dgncur_a(i,k,modefrm), 1.0e6*dgncur_a(i,k,modetoo)
 	    write(lunout,9820) 'xbeta ij0, ij3, ii0, jj0',   &
 		xbetaij0, xbetaij3, xbetaii0, xbetajj0
 	    write(lunout,9820) 'xbeta ij2i & j, ii2, jj2',   &
 		xbetaij2i, xbetaij2j, xbetaii2, xbetajj2
 	    write(lunout,9820) 'numbii, numbjj, deltat  ',   &
 		xnumbconc(modefrm), xnumbconc(modetoo), deltat
 	    write(lunout,9820) 'loss ij3, ii0, jj0      ',   &
 		(xbetaij3*xnumbconc(modetoo)*deltat),   &
 		(xbetaij0*xnumbconc(modetoo)*deltat+    &
 		 xbetaii0*xnumbconc(modefrm)*deltat),   &
 		(xbetajj0*xnumbconc(modetoo)*deltat)
 9801	format( / 72x, 'ACOAG' )
 9810	format( 'ACOAG ', a, 2i8, 3i7, 3(1pe15.6) )
 9820	format( 'ACOAG ', a, 4(1pe15.6) )
 9830	format( 'ACOAG ', a, i1, a, 4(1pe15.6) )
 	end if
 	end if
 	end if
 	end if   ! (ldiag2 > 0)
!   test diagnostics end ----------------------------------------------

	ybetaij0(ipair) = xbetaij0
	ybetaij3(ipair) = xbetaij3
	ybetaii0(ipair) = xbetaii0
	ybetajj0(ipair) = xbetajj0

	end do main_ipair1



	if ( (pair_option_acoag == 1) .or.   &
	     (pair_option_acoag == 2) ) then
!
!   calculate number and mass changes for pair_option_acoag == 1,2
!
main_ipair2: do ipair = 1, npair_acoag

	modefrm = modefrm_acoag(ipair)
	modetoo = modetoo_acoag(ipair)

!   calculate number changes
!   apply self-coagulation losses only once to a mode (when iselfcoagdone=0)
!	first calc change to "too" mode
!	next  calc change to "frm" mode, using average number conc of "too"
	if ( (mprognum_amode(modetoo) > 0) .and.   &
	     (iselfcoagdone(modetoo) <= 0) ) then
	    iselfcoagdone(modetoo) = 1
	    tmpn = xnumbconc(modetoo)
	    xnumbconcnew(modetoo) = tmpn/(1.0_r8 + deltat*ybetajj0(ipair)*tmpn)
	    xnumbconcavg(modetoo) = 0.5_r8*(xnumbconcnew(modetoo) + tmpn)
	    lstoo = numptr_amode(modetoo) - loffset
	    q(i,k,lstoo) = xnumbconcnew(modetoo)/aircon
	    dqdt(i,k,lstoo) = (xnumbconcnew(modetoo)-tmpn)*deltatinv_main/aircon
	end if

	if ( (mprognum_amode(modefrm) > 0) .and.   &
	     (iselfcoagdone(modefrm) <= 0) ) then
	    iselfcoagdone(modefrm) = 1
	    tmpn = xnumbconc(modefrm)
	    tmpa = deltat*ybetaij0(ipair)*xnumbconcavg(modetoo)
	    tmpb = deltat*ybetaii0(ipair)
	    tmpc = tmpa + tmpb*tmpn
	    if (abs(tmpc) < 0.01_r8) then
		xnumbconcnew(modefrm) = tmpn*exp(-tmpc)
	    else if (abs(tmpa) < 0.001_r8) then
		xnumbconcnew(modefrm) =   &
		    exp(-tmpa)*tmpn/(1.0_r8 + tmpb*tmpn)
	    else
		tmpf = tmpb*tmpn/tmpc
		tmpg = exp(-tmpa)
		tmph = tmpg*(1.0_r8 - tmpf)/(1.0_r8 - tmpg*tmpf)
		xnumbconcnew(modefrm) = tmpn*max( 0.0_r8, min( 1.0_r8, tmph ) )
	    end if
	    xnumbconcavg(modefrm) = 0.5_r8*(xnumbconcnew(modefrm) + tmpn)
	    lsfrm = numptr_amode(modefrm) - loffset
	    q(i,k,lsfrm) = xnumbconcnew(modefrm)/aircon
	    dqdt(i,k,lsfrm) = (xnumbconcnew(modefrm)-tmpn)*deltatinv_main/aircon
	end if

!   calculate mass changes
!     xbetaij3*xnumbconc(modetoo) = first order loss rate for modefrm volume
!     xferfracvol = fraction of modefrm volume transferred to modetoo over deltat
	dumloss = ybetaij3(ipair)*xnumbconcavg(modetoo)
	xferfracvol = 1.0_r8 - exp( -dumloss*deltat )
	xferfracvol = max( 0.0_r8, min( xferfrac_max, xferfracvol ) )

	do iq = 1, nspecfrm_acoag(ipair)
	    lsfrm = lspecfrm_acoag(iq,ipair) - loffset
	    lstoo = lspectoo_acoag(iq,ipair) - loffset
	    if (lsfrm > 0) then
		xferamt = q(i,k,lsfrm)*xferfracvol
		dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferamt*deltatinv_main
		q(i,k,lsfrm) = q(i,k,lsfrm) - xferamt
		if (lstoo > 0) then
		    dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferamt*deltatinv_main
		    q(i,k,lstoo) = q(i,k,lstoo) + xferamt
		end if
	    end if
	end do

	end do main_ipair2


	else if (pair_option_acoag == 3) then
!
!   calculate number and mass changes for pair_option_acoag == 3
!

!   calculate number changes to accum mode
	if (mprognum_amode(macc) > 0) then
	    tmpn = xnumbconc(macc)
	    xnumbconcnew(macc) = tmpn/(1.0_r8 + deltat*ybetajj0(ip_aitacc)*tmpn)
	    xnumbconcavg(macc) = 0.5_r8*(xnumbconcnew(macc) + tmpn)
	    lstoo = numptr_amode(macc) - loffset
	    q(i,k,lstoo) = xnumbconcnew(macc)/aircon
	    dqdt(i,k,lstoo) = (xnumbconcnew(macc)-tmpn)*deltatinv_main/aircon
	end if

!   calculate number changes to primary carbon mode
	modefrm = modeptr_pcarbon
	if (mprognum_amode(mpca) > 0) then
	    tmpn = xnumbconc(mpca)
	    tmpa = deltat*ybetaij0(ip_pcaacc)*xnumbconcavg(macc)
	    tmpb = deltat*ybetaii0(ip_pcaacc)
	    tmpc = tmpa + tmpb*tmpn
	    if (abs(tmpc) < 0.01_r8) then
		xnumbconcnew(mpca) = tmpn*exp(-tmpc)
	    else if (abs(tmpa) < 0.001_r8) then
		xnumbconcnew(mpca) =   &
		    exp(-tmpa)*tmpn/(1.0_r8 + tmpb*tmpn)
	    else
		tmpf = tmpb*tmpn/tmpc
		tmpg = exp(-tmpa)
		tmph = tmpg*(1.0_r8 - tmpf)/(1.0_r8 - tmpg*tmpf)
		xnumbconcnew(mpca) = tmpn*max( 0.0_r8, min( 1.0_r8, tmph ) )
	    end if
	    xnumbconcavg(mpca) = 0.5_r8*(xnumbconcnew(mpca) + tmpn)
	    lsfrm = numptr_amode(mpca) - loffset
	    q(i,k,lsfrm) = xnumbconcnew(mpca)/aircon
	    dqdt(i,k,lsfrm) = (xnumbconcnew(mpca)-tmpn)*deltatinv_main/aircon
	end if

!   calculate number changes to aitken mode
	if (mprognum_amode(mait) > 0) then
	    tmpn = xnumbconc(mait)
	    tmpa = deltat*( ybetaij0(ip_aitacc)*xnumbconcavg(macc)   &
	                  + ybetaij0(ip_aitpca)*xnumbconcavg(mpca) )
	    tmpb = deltat*ybetaii0(ip_aitacc)
	    tmpc = tmpa + tmpb*tmpn
	    if (abs(tmpc) < 0.01_r8) then
		xnumbconcnew(mait) = tmpn*exp(-tmpc)
	    else if (abs(tmpa) < 0.001_r8) then
		xnumbconcnew(mait) =   &
		    exp(-tmpa)*tmpn/(1.0_r8 + tmpb*tmpn)
	    else
		tmpf = tmpb*tmpn/tmpc
		tmpg = exp(-tmpa)
		tmph = tmpg*(1.0_r8 - tmpf)/(1.0_r8 - tmpg*tmpf)
		xnumbconcnew(mait) = tmpn*max( 0.0_r8, min( 1.0_r8, tmph ) )
	    end if
	    xnumbconcavg(mait) = 0.5_r8*(xnumbconcnew(mait) + tmpn)
	    lsfrm = numptr_amode(mait) - loffset
	    q(i,k,lsfrm) = xnumbconcnew(mait)/aircon
	    dqdt(i,k,lsfrm) = (xnumbconcnew(mait)-tmpn)*deltatinv_main/aircon
	end if


!   calculate mass changes from aitken-->accum direct coagulation and
!	aitken-->pcarbon-->accum coagulation/aging
!   also calc volume of shell material (so4 & nh4 from aitken-->pcarbon)
	dumloss = ybetaij3(ip_aitacc)*xnumbconcavg(macc)   &
	        + ybetaij3(ip_aitpca)*xnumbconcavg(mpca)
	tmpa = ybetaij3(ip_aitpca)*xnumbconcavg(mpca)/max( dumloss, 1.0e-37_r8 )
	xferfracvol = 1.0_r8 - exp( -dumloss*deltat )
	xferfracvol = max( 0.0_r8, min( xferfrac_max, xferfracvol ) )
	vol_shell = 0.0

	ipair = ip_aitacc
	do iq = 1, nspecfrm_acoag(ipair)
	    lsfrm = lspecfrm_acoag(iq,ipair) - loffset
	    lstoo = lspectoo_acoag(iq,ipair) - loffset
	    if (lsfrm > 0) then
		xferamt = q(i,k,lsfrm)*xferfracvol
		dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferamt*deltatinv_main
		q(i,k,lsfrm) = q(i,k,lsfrm) - xferamt
		if (lstoo > 0) then
		    dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferamt*deltatinv_main
		    q(i,k,lstoo) = q(i,k,lstoo) + xferamt
		end if
		vol_shell = vol_shell + xferamt*tmpa*fac_m2v_aitage(iq)
	    end if
	end do


!   now calculate aging transfer fraction for pcarbon-->accum
!   this duplicates the code in modal_aero_gasaerexch
	vol_core = 0.0
	do l = 1, nspec_amode(mpca)
	    vol_core = vol_core + &
		q(i,k,lmassptr_amode(l,mpca)-loffset)*fac_m2v_pcarbon(l)
	end do
	tmp1 = vol_shell*dgncur_a(i,k,mpca)*fac_volsfc_pcarbon
	tmp2 = 6.0_r8*dr_so4_monolayers_pcage*vol_core
	tmp2 = max( tmp2, 0.0_r8 )
	if (tmp1 >= tmp2) then
	    xferfrac_pcage = xferfrac_max
	else
	    xferfrac_pcage = min( tmp1/tmp2, xferfrac_max )
	end if


!   calculate mass changes from pcarbon-->accum by direct coagulation
!   and aging
	dumloss = ybetaij3(ip_pcaacc)*xnumbconcavg(macc)
	xferfracvol = 1.0_r8 - exp( -dumloss*deltat )
	xferfracvol = xferfracvol + xferfrac_pcage
	xferfracvol = max( 0.0_r8, min( xferfrac_max, xferfracvol ) )

	ipair = ip_pcaacc
	do iq = 1, nspecfrm_acoag(ipair)
	    lsfrm = lspecfrm_acoag(iq,ipair) - loffset
	    lstoo = lspectoo_acoag(iq,ipair) - loffset
	    if (lsfrm > 0) then
		xferamt = q(i,k,lsfrm)*xferfracvol
		dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferamt*deltatinv_main
		q(i,k,lsfrm) = q(i,k,lsfrm) - xferamt
		if (lstoo > 0) then
		    dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferamt*deltatinv_main
		    q(i,k,lstoo) = q(i,k,lstoo) + xferamt
		end if
	    end if
	end do

	lsfrm = numptr_amode(mpca) - loffset
	lstoo = numptr_amode(macc) - loffset
	if (lsfrm > 0) then
	    xferamt = q(i,k,lsfrm)*xferfrac_pcage
	    dqdt(i,k,lsfrm) = dqdt(i,k,lsfrm) - xferamt*deltatinv_main
	    q(i,k,lsfrm) = q(i,k,lsfrm) - xferamt
	    if (lstoo > 0) then
		dqdt(i,k,lstoo) = dqdt(i,k,lstoo) + xferamt*deltatinv_main
		q(i,k,lstoo) = q(i,k,lstoo) + xferamt
	    end if
	end if



	else   ! (pair_option_acoag /= 1,2,3) then

	write(lunout,*) '*** modal_aero_coag_sub error'
	write(lunout,*) '    cannot do _coag_sub error pair_option_acoag =', &
		pair_option_acoag
	call endrun( 'modal_aero_coag_sub error' )


	end if   ! (pair_option_acoag == ...)


!   test diagnostics begin --------------------------------------------
 	if (ldiag3 > 0) then
 	if (nstep <= 3) then
 	if ((lonndx(i) == 37) .and. (latndx(i) == 23)) then
 	if ((mod(k-1,5) == 0) .or. (k>=23)) then
 	   if (pair_option_acoag == 3) then
 		write(*,*)
 		write(lunout,9820) 'xnumbconcavg ait,acc,pca', &
 		    xnumbconcavg(mait), xnumbconcavg(macc), xnumbconcavg(mpca)
 		write(lunout,9820) 'vshell, core            ', &
 		    vol_shell, vol_core
 		write(lunout,9820) 'dr_mono, dgn            ', &
 		    dr_so4_monolayers_pcage, dgncur_a(i,k,mpca)
 		write(lunout,9820) 'tmp1, tmp2              ', tmp1, tmp2
 		write(lunout,9820) 'xferfrac_age            ', xferfrac_pcage
 	   end if

 	   do ipair = 1, npair_acoag
 	   modefrm = modefrm_acoag(ipair)
 	   modetoo = modetoo_acoag(ipair)
 	   if (npair_acoag > 1) then
 		write(lunout,*)
 		write(lunout,9810) 'ipair =   ', ipair
 	   end if

 	   do iq = 1, nspecfrm_acoag(ipair)
 	   lsfrm = lspecfrm_acoag(iq,ipair) - loffset
 	   lstoo = lspectoo_acoag(iq,ipair) - loffset
 	   if (lsfrm > 0) then
 	   tmp_qold = q(i,k,lsfrm) - dqdt(i,k,lsfrm)*deltat_main
!	   write(lunout,9820) 'm1 frm dqdt/q0,dqdt,q0/1',   &
 	   write(lunout,9830) 'm', iq,   &
 	                        ' frm dqdt/q0,dqdt,q0/1',   &
 		dqdt(i,k,lsfrm)/tmp_qold, dqdt(i,k,lsfrm), tmp_qold, q(i,k,lsfrm)
 	   end if
 	   if (lstoo > 0) then
 	   tmp_qold = q(i,k,lstoo) - dqdt(i,k,lstoo)*deltat_main
 	   write(lunout,9830) 'm', iq,   &
 	                        ' too dqdt/q0,dqdt,q0/1',   &
 		dqdt(i,k,lstoo)/tmp_qold, dqdt(i,k,lstoo), tmp_qold, q(i,k,lstoo)
 	   end if
 	   end do   ! iq

 	   lsfrm = numptr_amode(modefrm) - loffset
 	   lstoo = numptr_amode(modetoo) - loffset
 	   if (lsfrm > 0) then
 	   tmp_qold = q(i,k,lsfrm) - dqdt(i,k,lsfrm)*deltat_main
 	   write(lunout,9820) 'n  frm dqdt/q0,dqdt,q0/1',   &
 		dqdt(i,k,lsfrm)/tmp_qold, dqdt(i,k,lsfrm), tmp_qold, q(i,k,lsfrm)
 	   end if
 	   if (lstoo > 0) then
 	   tmp_qold = q(i,k,lstoo) - dqdt(i,k,lstoo)*deltat_main
 	   write(lunout,9820) 'n  too dqdt/q0,dqdt,q0/1',   &
 		dqdt(i,k,lstoo)/tmp_qold, dqdt(i,k,lstoo), tmp_qold, q(i,k,lstoo)
 	   end if

 	   end do   ! ipair
 	end if
 	end if
 	end if
 	end if   ! (ldiag3 > 0)
!   test diagnostics end ----------------------------------------------



	end do main_i
	end do main_k


! set dotend's
	do ipair = 1, npair_acoag
	    modefrm = modefrm_acoag(ipair)
	    modetoo = modetoo_acoag(ipair)

	    do iq = 1, nspecfrm_acoag(ipair)
		lsfrm = lspecfrm_acoag(iq,ipair) - loffset
		lstoo = lspectoo_acoag(iq,ipair) - loffset
		if (lsfrm > 0) dotend(lsfrm) = .true.
		if (lstoo > 0) dotend(lstoo) = .true.
	    end do

	    if (mprognum_amode(modefrm) > 0) then
		lsfrm = numptr_amode(modefrm) - loffset
		if (lsfrm > 0) dotend(lsfrm) = .true.
	    end if
	    if (mprognum_amode(modetoo) > 0) then
		lstoo = numptr_amode(modetoo) - loffset
		if (lstoo > 0) dotend(lstoo) = .true.
	    end if

	end do


!   do history file column-tendency fields
	do l = loffset+1, pcnst
	    lmz = l - loffset
	    if ( .not. dotend(lmz) ) cycle

	    qsrflx(:) = 0.0_r8
	    do k = 1, pver
	    do i = 1, ncol
		qsrflx(i) = qsrflx(i) + dqdt(i,k,lmz)*pdel(i,k)
	    end do
	    end do
	    qsrflx(:) = qsrflx(:)*(adv_mass(lmz)/(gravit*mwdry))
	    fieldname = trim(cnst_name(l)) // '_sfcoag1'
	    call outfld( fieldname, qsrflx, pcols, lchnk )
!	    if (( masterproc ) .and. (nstep < 1)) &
!		write(*,'(2(a,2x),1p,e11.3)') &
!		'modal_aero_coag_sub outfld', fieldname, adv_mass(lmz)
	end do ! l = ...


	return


!EOC
	end subroutine modal_aero_coag_sub


!----------------------------------------------------------------------
!----------------------------------------------------------------------
	subroutine modal_aero_coag_init
!
!   computes pointers for species transfer during coagulation
!
	use modal_aero_data
	use modal_aero_gasaerexch, only:  &
		modefrm_pcage, nspecfrm_pcage, lspecfrm_pcage, lspectoo_pcage

	use abortutils,      only: endrun
	use cam_history,     only: addfld, add_default, fieldname_len, phys_decomp
	use constituents,    only: pcnst, cnst_name
	use spmd_utils,      only: masterproc
        use phys_control,    only: phys_getopts

	implicit none

!   local variables
	integer :: ipair, iq, iqfrm, iqfrm_aa, iqtoo, iqtoo_aa
	integer :: l, lsfrm, lstoo, lunout
	integer :: m, mfrm, mtoo, mtef
	integer :: nsamefrm, nsametoo, nspec

	character(len=fieldname_len)   :: tmpname
	character(len=fieldname_len+3) :: fieldname
	character(128)                 :: long_name
	character(8)                   :: unit

	logical :: dotend(pcnst)
        logical :: history_aerosol      ! Output the MAM aerosol tendencies
 
        !-----------------------------------------------------------------------     
    
        call phys_getopts( history_aerosol_out        = history_aerosol   )

	lunout = 6
!
!   define "from mode" and "to mode" for each coagulation pairing
!	currently just a2-->a1 coagulation
!
	if (pair_option_acoag == 1) then
	    npair_acoag = 1
	    modefrm_acoag(1) = modeptr_aitken
	    modetoo_acoag(1) = modeptr_accum
	    modetooeff_acoag(1) = modeptr_accum
	else if (pair_option_acoag == 2) then
	    npair_acoag = 2
	    modefrm_acoag(1) = modeptr_aitken
	    modetoo_acoag(1) = modeptr_accum
	    modetooeff_acoag(1) = modeptr_accum
	    modefrm_acoag(2) = modeptr_pcarbon
	    modetoo_acoag(2) = modeptr_accum
	    modetooeff_acoag(2) = modeptr_accum
	else if (pair_option_acoag == 3) then
	    npair_acoag = 3
	    modefrm_acoag(1) = modeptr_aitken
	    modetoo_acoag(1) = modeptr_accum
	    modetooeff_acoag(1) = modeptr_accum
	    modefrm_acoag(2) = modeptr_pcarbon
	    modetoo_acoag(2) = modeptr_accum
	    modetooeff_acoag(2) = modeptr_accum
	    modefrm_acoag(3) = modeptr_aitken
	    modetoo_acoag(3) = modeptr_pcarbon
	    modetooeff_acoag(3) = modeptr_accum
	    if (modefrm_pcage <= 0) then
		write(*,*) '*** modal_aero_coag_init error'
		write(*,*) '    pair_option_acoag, modefrm_pcage mismatch'
		write(*,*) '    pair_option_acoag, modefrm_pcage =', &
		    pair_option_acoag, modefrm_pcage
		call endrun( 'modal_aero_coag_init error' )
	    end if
	else
	    npair_acoag = 0
	    return
	end if

!
!   define species involved in each coagulation pairing
!	(include aerosol water)
!
aa_ipair: do ipair = 1, npair_acoag

	mfrm = modefrm_acoag(ipair)
	mtoo = modetoo_acoag(ipair)
	mtef = modetooeff_acoag(ipair)
	if ( (mfrm < 1) .or. (mfrm > ntot_amode) .or.   &
	     (mtoo < 1) .or. (mtoo > ntot_amode) .or.   &
	     (mtef < 1) .or. (mtef > ntot_amode) ) then
	    write(*,*) '*** modal_aero_coag_init error'
	    write(*,*) '    ipair, ntot_amode =', ipair, ntot_amode
	    write(*,*) '    mfrm, mtoo, mtef  =', mfrm, mtoo, mtef
	    call endrun( 'modal_aero_coag_init error' )
	end if


	mtoo = mtef   ! effective modetoo
	nspec = 0
aa_iqfrm: do iqfrm = 1, nspec_amode(mfrm)
	    lsfrm = lmassptr_amode(iqfrm,mfrm)
	    if ((lsfrm .lt. 1) .or. (lsfrm .gt. pcnst)) cycle aa_iqfrm

! find "too" species having same lspectype_amode as the "frm" species
! several species in a mode may have the same lspectype_amode, so also
!    use the ordering as a criterion (e.g., 1st <--> 1st, 2nd <--> 2nd)
	    iqfrm_aa = 1
	    iqtoo_aa = 1
	    if (iqfrm .gt. nspec_amode(mfrm)) then
		iqfrm_aa = nspec_amode(mfrm) + 1
		iqtoo_aa = nspec_amode(mtoo) + 1
	    end if
	    nsamefrm = 0
	    do iq = iqfrm_aa, iqfrm
		if ( lspectype_amode(iq   ,mfrm) .eq.   &
      		     lspectype_amode(iqfrm,mfrm) ) then
		    nsamefrm = nsamefrm + 1
		end if
	    end do
	    nsametoo = 0
	    lstoo = 0
	    do iqtoo = iqtoo_aa, nspec_amode(mtoo)
		if ( lspectype_amode(iqtoo,mtoo) .eq.   &
      		     lspectype_amode(iqfrm,mfrm) ) then
		    nsametoo = nsametoo + 1
		    if (nsametoo .eq. nsamefrm) then
			lstoo = lmassptr_amode(iqtoo,mtoo)
			exit
		    end if
		end if
	    end do

	    nspec = nspec + 1
	    lspecfrm_acoag(nspec,ipair) = lsfrm
	    lspectoo_acoag(nspec,ipair) = lstoo
	end do aa_iqfrm

!       lsfrm = lwaterptr_amode(mfrm)
!       if ((lsfrm .ge. 1) .and. (lsfrm .le. pcnst)) then
!           lstoo = lwaterptr_amode(mtoo)
!           if ((lstoo .lt. 1) .or. (lstoo .gt. pcnst)) lstoo = 0
!           nspec = nspec + 1
!           lspecfrm_acoag(nspec,ipair) = lsfrm
!           lspectoo_acoag(nspec,ipair) = lstoo
!       end if

	nspecfrm_acoag(ipair) = nspec
	end do aa_ipair

!
!   output results
!
	if ( masterproc ) then

	write(lunout,9310)

	do ipair = 1, npair_acoag
	  mfrm = modefrm_acoag(ipair)
	  mtoo = modetoo_acoag(ipair)
	  mtef = modetooeff_acoag(ipair)
	  write(lunout,9320) ipair, mfrm, mtoo, mtef

	  do iq = 1, nspecfrm_acoag(ipair)
	    lsfrm = lspecfrm_acoag(iq,ipair)
	    lstoo = lspectoo_acoag(iq,ipair)
	    if (lstoo .gt. 0) then
		write(lunout,9330) lsfrm, cnst_name(lsfrm),   &
      			lstoo, cnst_name(lstoo)
	    else
		write(lunout,9340) lsfrm, cnst_name(lsfrm)
	    end if
	  end do

	end do ! ipair = ...
	write(lunout,*)

	end if ! ( masterproc ) 

9310	format( / 'subr. modal_aero_coag_init' )
9320	format( 'pair', i3, 5x, 'mode', i3, &
		' ---> mode', i3, '   eff', i3 )
9330	format( 5x, 'spec', i3, '=', a, ' ---> spec', i3, '=', a )
9340	format( 5x, 'spec', i3, '=', a, ' ---> LOSS' )


!
!   create history file column-tendency fields
!
	dotend(:) = .false.
	do ipair = 1, npair_acoag
	  do iq = 1, nspecfrm_acoag(ipair)
	    l = lspecfrm_acoag(iq,ipair)
	    if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
	    l = lspectoo_acoag(iq,ipair)
	    if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
	  end do

	  m = modefrm_acoag(ipair)
	  if ((m > 0) .and. (m <= ntot_amode)) then
	    l = numptr_amode(m)
	    if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
	  end if
	  m = modetoo_acoag(ipair)
	  if ((m > 0) .and. (m <= ntot_amode)) then
	    l = numptr_amode(m)
	    if ((l > 0) .and. (l <= pcnst)) dotend(l) = .true.
	  end if
	end do ! ipair = ...

	if (pair_option_acoag == 3) then
	   do iq = 1, nspecfrm_pcage
	      lsfrm = lspecfrm_pcage(iq)
	      lstoo = lspectoo_pcage(iq)
	      if ((lsfrm > 0) .and. (lsfrm <= pcnst)) then
	         dotend(lsfrm) = .true.
	         if ((lstoo > 0) .and. (lstoo <= pcnst)) then
	            dotend(lstoo) = .true.
	         end if
	      end if
	   end do
	end if

	do l = 1, pcnst
	    if ( .not. dotend(l) ) cycle
	    tmpname = cnst_name(l)
	    unit = 'kg/m2/s'
	    do m = 1, ntot_amode
	        if (l == numptr_amode(m)) unit = '#/m2/s'
	    end do
	    fieldname = trim(tmpname) // '_sfcoag1'
	    long_name = trim(tmpname) // ' modal_aero coagulation column tendency'
	    call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
            if ( history_aerosol ) then 
               call add_default( fieldname, 1, ' ' )
	    endif
	    if ( masterproc ) write(*,'(3(a,2x))') &
		'modal_aero_coag_init addfld', fieldname, unit
	end do ! l = ...


	return
	end subroutine modal_aero_coag_init


!----------------------------------------------------------------------
!----------------------------------------------------------------------
	subroutine calc_coag_coef4(   &
          dgni, dgnj, alnsgi, alnsgj, rhopi, rhopj,   &
          wetddrydi, wetddrydj,   &
          temp, presscgs, lunerr, lunout, iok,   &
          betaij0, betaij2i, betaij2j, betaij3,   &
          betaii0, betaii2, betajj0, betajj2 )
!
!   computes the following coagulation rate "coefficients"
!   for lognormal modes
!
!   self coagulation
!	dNi/dt = - betaii0*Ni*Ni
!	dNj/dt = - betajj0*Nj*Nj
!	dSi/dt = - betaii2*Si*Ni
!	dSj/dt = - betajj2*Sj*Nj
!
!   modei-modej coagulation
!	dNi/dt = - betaij0*Ni*Nj
!	dNj/dt = 0.
!	dSi/dt = - betaij2i*Si*Nj
!	dSj/dt = + betaij2j*Si*Nj
!	dVi/dt  = - betaij3*Vi*Nj == - dVj/dt
!
!   Ni, Nj are number  concentrations (particles/cm3)
!   Si, Sj are surface concentrations (cm2/cm3)
!   Vi, Vj are volume  concentrations (cm3/cm3)
!   rates are (N,S,V)/s
!
!   mode i is the smaller mode (e.g., Aitken)
!   mode j is the larger  mode (e.g., accumulation)
!
!   input arguments
!	dgni, dgnj = DRY median diameters for number distribution (cm)
!	alnsgi, alnsgj = ln of geometric standard deviations (dimensionless)
!	rhopi, rhopj = WET density (g/cm3)
!	wetddrydi, wetddrydj = (WET median diameter)/(DRY median diameter)
!	temp = air temperature (K)
!	presscgs = air pressure (dynes/cm2)
!	lunerr, lunout = logical unit for error or diagnostic output
!
!   output arguments
!	iok = status flag  (+1 = success, 0/negative = failure)
!	beta--- = see above
!

	implicit none

!   arguments
	integer, intent(in)  :: lunerr, lunout
	integer, intent(out) :: iok
	real(r4), intent(in) ::   &
          dgni, dgnj, alnsgi, alnsgj, rhopi, rhopj,   &
          wetddrydi, wetddrydj,   &
          temp, presscgs
	real(r4), intent(out) ::   &
          betaij0, betaij2i, betaij2j, betaij3,   &
          betaii0, betaii2, betajj0, betajj2

!   local variables
	real(r4) airprs, airtemp
	real(r4) dgacc, dgatk, pdensac, pdensat
	real(r4) batat(2), batac(2), bacat(2), bacac(2), c3ij
	real(r4) dp2bar_mks_i, dp2bar_mks_j, dp3bar_mks_i

!   check for reasonable inputs
	iok = -1
	if ((dgni .lt. 1.e-8) .or. (dgni .gt. 1.)) then
	    write(lunerr,9100) 'dgni', dgni
	    return
	else if ((dgnj .lt. 1.e-8) .or. (dgnj .gt. 1.)) then
	    write(lunerr,9100) 'dgnj', dgnj
	    return
!	else if (dgni .gt. dgnj) then
!	    write(lunerr,9110) dgni, dgnj
!	    return
	else if ((alnsgi .lt. 0.0) .or. (alnsgi .gt. 2.3)) then
	    write(lunerr,9100) 'alnsgi', alnsgi
	    return
	else if ((alnsgj .lt. 0.0) .or. (alnsgj .gt. 2.3)) then
	    write(lunerr,9100) 'alnsgj', alnsgj
	    return
	else if ((rhopi .lt. 0.01) .or. (rhopi .gt. 100.)) then
	    write(lunerr,9100) 'rhopi', rhopi
	    return
	else if ((rhopj .lt. 0.01) .or. (rhopj .gt. 100.)) then
	    write(lunerr,9100) 'rhopj', rhopj
	    return
	else if ((temp .lt. 10.) .or. (temp .gt. 370.)) then
	    write(lunerr,9100) 'temp', temp
	    return
	else if ((presscgs .lt. 1.e2) .or. (presscgs .gt. 1.5e6)) then
	    write(lunerr,9100) 'presscgs', presscgs
	    return
	end if
9100	format( '*** subr. calc_coag_coef_4 - bad value for ', a,   &
      		' = ', 1pe15.5 )
!9110	format( '*** subr. calc_coag_coef_4 - dgni > dgnj -- ',
!     +		2(1pe15.5) )
	iok = +1

!
!   define mks variables
!
	airtemp = temp
!   dyne/cm2 to pa
	airprs = presscgs * 1.0e-1

!   cm to m
	dgatk = dgni * 1.0e-2
	dgacc = dgnj * 1.0e-2

!   g/cm3 to kg/m3
	pdensat = rhopi * 1.0e+3
	pdensac = rhopj * 1.0e+3

!
!   call interace to binkowski routines
!
	call bink_coag_rates(   &
          dgatk, dgacc, alnsgi, alnsgj, pdensat, pdensac,   &
          wetddrydi, wetddrydj,   &
          airtemp, airprs, lunerr, iok,   &
          batat, batac, bacat, bacac, c3ij )

!
!   transfer the batat, batac, bacat, bacac values
!   to the beta--- variables
!

!   self coagulation, number
!	dNi/dt = - betaii0*Ni*Ni
!	dNXi/dt = - batat(1)*NXi*NXi
!   where Ni is (particles/cm3), NXi is (particles/m3)
!   the first-order loss rates are just (1/s) and thus are equal
!	d[ln(Ni)]/dt = d[ln(NXi)]/dt
!   so
!	betaii0 = batat(1) * (NXi/Ni)
!   conversion factors are
!	1.0e-20 to undo 1.0e+20 scaling done in intercoag & intracoag
!	(NXi/Ni) = 1.0e6
!
	betaii0 = batat(1) * 1.0e-14
	betajj0 = bacac(1) * 1.0e-14

!   self coagulation, surface
!	dSi/dt = - betaii2*Si*Ni
!	dM2Xi/dt = - batat(1)*NXi*NXi
!   where Si is surface in (cm2/cm3) and M2Xi is (surface/pi) in (m2/m3)
!   the first-order loss rates are just (1/s) and thus are equal
!	betaii2 = batat(2) * (NXi/M2Xi) * (NXi/Ni)
!   conversion factors are
!	1.0e-20 to undo 1.0e+20 scaling done in intercoag & intracoag
!	(NXi/Ni) = 1.0e6
!	(M2Xi/NXi) = dp2bar_mks_i/pi, where dp2bar_mks_i is the average
!	    Dp**2 in m**2
!
	dp2bar_mks_j = (dgnj**2) * exp( 2.0*alnsgj*alnsgj ) * 1.0e-4

	dp2bar_mks_i = (dgni**2) * exp( 2.0*alnsgi*alnsgi ) * 1.0e-4
	dp3bar_mks_i = (dgni**3) * exp( 4.5*alnsgi*alnsgi ) * 1.0e-6

	betaii2 = (batat(2) / dp2bar_mks_i) * 1.0e-14
	betajj2 = (bacac(2) / dp2bar_mks_j) * 1.0e-14

!   modei-modej coagulation, number
!	dNi/dt = - betaij0*Ni*Nj
!	dNj/dt = 0.
!   conversions as for self coagulation
!
	betaij0  = batac(1) * 1.0e-14

!   modei-modej coagulation, surface
!	dSi/dt = - betaij2i*Si*Nj
!	dSj/dt = + betaij2j*Si*Nj
!   conversions as for self coagulation
!
	betaij2i = (batac(2) / dp2bar_mks_i) * 1.0e-14
	betaij2j = (bacat(2) / dp2bar_mks_i) * 1.0e-14

!   modei-modej coagulation, volume
!	dVi/dt  = - betaij3*Vi*Nj
!	dM3Xi/dt = - c3ij*NXi*NXj
!   where Vi is (cm3/cm3) and M3Xi is (volume*6/pi) in (m3/m3)
!   the first-order loss rates are just (1/s) and thus are equal
!	betaij3 = c3ij * (NXi/M3Xi) * (NXj/Nj)
!   conversion factors are
!	1.0e-20 to undo 1.0e+20 scaling done in intercoag & intracoag
!	(NXj/Nj) = 1.0e6
!	(M3Xi/NXi) = dp3bar_mks_i*(6/pi), where dp3bar_mks_i is the average
!	    Dp**3 in m**3
!
	betaij3  = ( c3ij / dp3bar_mks_i ) * 1.0e-14

	return
	end subroutine calc_coag_coef4



!-----------------------------------------------------------------------
	subroutine bink_coag_rates(   &
          dgatk, dgacc, xxlsgat, xxlsgac, pdensat, pdensac,   &
          wetddrydat, wetddrydac,   &
          airtemp, airprs, lunerr, iok,   &
          batat, batac, bacat, bacac, c3ij )

!
!   provides interface to F. Binkowski's intracoag and intercoag
!
!   this code was "cut" from F. Binkowski's aero_info_ae3.f & aero_subs3_ae3.f
!
!   computes the following coagulation rate "coefficients"
!
!   self coagulation
!	dNXi/dt  = - (1.e-20*batat(1))*NXi*NXi
!	dNXj/dt  = - (1.e-20*bacac(1))*NXj*NXj
!	dM2Xi/dt = - (1.e-20*batat(2))*NXi*NXi
!	dM2Xj/dt = - (1.e-20*bacac(2))*NXj*NXj
!
!   modei-modej coagulation
!	dNXi/dt  = - (1.e-20*batac(1))*NXi*NXj
!	dNXj/dt  = 0.
!	dM2Xi/dt = - (1.e-20*batac(2))*NXi*NXj
!	dM2Xj/dt = + (1.e-20*bacat(2))*NXi*NXj
!	dM3Xi/dt =     - (1.e-20*c3ij)*NXi*NXj == - dM3Xj/dt
!
!   NXi,  NXj are number concentrations (particles/m3)
!   M2Xi, M2Xj are 2nd moment (surface/pi)  concentrations (m2/m3)
!   M3Xi, M3Xj are 2rd moment (volume*6/pi) concentrations (m3/m3)
!
!   in the above, mode i is aitken mode, mode j is accumulation mode
!
!   *** note that the batat, batac, bacat, bacac, c3ij
!	all must be multiplied by 1.e-20 to undo the 1.e+20 scaling
!	in the real*4 versions of intracoag and intercoag
!
!   input arguments
!	dgatk, dgacc = DRY median diameters for number distribution (m)
!	    for aitken (atk or at) and accumulation (acc or ac) modes
!	xxlsgat, xxlsgac = ln of geometric standard deviations (dimensionless)
!	pdensat, pdensac = WET density (kg/m3)
!	wetddrydat, wetddrydac = (WET median diameter)/(DRY median diameter)
!	airtemp = air temperature (K)
!	airprs = air pressure (Pa)
!	lunerr  = logical unit for error output
!
!   output arguments
!	iok = status flag  (+1 = success, 0/negative = failure)
!	batat, batac, bacat, bacac, c3ij = see above
!
	implicit none

!   arguments
	integer lunerr, iok
	real(r4) dgatk, dgacc, xxlsgat, xxlsgac, pdensat, pdensac,   &
          wetddrydat, wetddrydac,   &
          airtemp, airprs
	real(r4) batat(2), batac(2), bacat(2), bacac(2), c3ij

! *** modal geometric mean diameters: [ m ]
!     real dgatk   ! nuclei mode
!     real dgacc   ! accumulation mode
! *** log of modal geometric standard deviation
!     real xxlsgat  ! aitken mode
!     real xxlsgac  ! accumulation mode
! *** average modal particle densities  [ kg/m**3 ]
!     real pdensat         ! nuclei mode
!     real pdensac         ! accumulation mode

!     real airtemp   ! air temperature [ k ]
!     real airprs    ! air pressure in [ pa ]


!   local variables

      real(r4)   two3
      parameter( two3 = 2.0/3.0 )

      real(r4)    avo      ! avogadro's constant [ 1/mol ]
      parameter ( avo = 6.0221367e23 )

      real(r4)    rgasuniv ! universal gas constant [ j/mol-k ]
      parameter ( rgasuniv = 8.314510 )

      real(r4)    boltz    ! boltzmann's constant [ j / k]
      parameter ( boltz = rgasuniv / avo )

      real(r4)    p0       !  starting standard surface pressure [ pa ]
      parameter ( p0 = 101325.0 )

      real(r4)    t0       !  starting standard surface temperature [ k ]
      parameter ( t0 = 288.15 )

      real(r4) xlm    ! atmospheric mean free path [ m ]
      real(r4) amu    ! atmospheric dynamic viscosity [ kg/m s ]

      real(r4) kfm, knc, lamda, sqrt_temp


! fsb calculate the square root of the ambient
!      temperature for later use

! *** calculate mean free path [ m ]:
! *** 6.6328e-8 is the sea level values given in table i.2.8
! *** on page 10 of u.s. standard atmosphere 1962
      xlm = 6.6328e-8 * p0 * airtemp / ( t0 * airprs )

! ***   calculate dynamic viscosity [ kg m**-1 s**-1 ]:
! *** u.s. standard atmosphere 1962 page 14 expression
!     for dynamic viscosity is:
!     dynamic viscosity =  beta * t * sqrt(t) / ( t + s)
!     where beta = 1.458e-6 [ kg sec^-1 k**-0.5 ], s = 110.4 [ k ].
      sqrt_temp = sqrt( airtemp)
      amu = 1.458e-6 * airtemp * sqrt_temp / ( airtemp + 110.4 )



! *** coagulation
! *** set up coagulation rates

! *** moment independent factors
            knc = two3 * boltz *  airtemp / amu
            lamda = xlm

! ***  calculate the coagulation coefficients for use
!      in the  aitken (nuclei) & accumulation modes
! ***  with gauss-hermite numerical quadrature
!      using 10 abscissas.

! *** aitken - aitken mode coagulation
      kfm = sqrt( 3.0 * boltz * airtemp / pdensat )
      call  intracoag_gh(lamda,   &
                         kfm, knc,   &
                         dgatk,   &
                         xxlsgat,   &
                         wetddrydat,   &
                         batat(2),   &
                         batat(1) )

! *** accumulation - accumulation mode coagulation
      kfm = sqrt( 3.0 * boltz * airtemp / pdensac )
      call  intracoag_gh(lamda,   &
                         kfm, knc,   &
                         dgacc,   &
                         xxlsgac,   &
                         wetddrydac,   &
                         bacac(2),   &
                         bacac(1) )

! *** aitken accumulation mode coagulation
      bacat(1) = 0.0 ! not used
      kfm  = sqrt( 6.0 * boltz * airtemp  /   &
                             ( pdensat + pdensac ) )
      call  intercoag_gh(lamda,   &
                         kfm, knc,   &
                         dgatk, dgacc ,   &
                         xxlsgat, xxlsgac ,   &
                         wetddrydat, wetddrydac,   &
                         batac(2),   &
                         bacat(2),   &
                         batac(1),   &
                         c3ij          )

!     c30atac = c3ij * cblk( vac0 ) * cblk( vat0 )
!     loss =  c30atac / cblk( vat3 )

      return
      end subroutine bink_coag_rates



! ---------------------------------------------------------------------
! fsb subrs to do gauss-hermite numerical quadrature

      subroutine intracoag_gh( lamda, kfm, knc,   &
           dg, xlnsig, wetddryd,   &
           quads11, quadn11)

!  fsb this version is for intramodal coagulation for number
!       and second moment

! fsb this version runs in real*4 arithmetic

! *** this version calculates the coagulation coefficients
!     using the harmonic mean approach for both fm and nc cases.
! *** does gauss-hermite quadrature for intra-modal
!      coagulation integrals for 2nd moment
!     for a lognormal distribution
!      defined by dg,xlnsig,
! *** dg and xlnsig are the geometric mean diameters (meters)
!      and the logarithms of the
!     geometric standard deviations (dimensionless)
!     whose meaning is defined below at the end of the routine
!     ghxi, ghwi are the gauss-hermite weights and n is one-half the
!     number of abscissas, since an even number of abscissas is used
!
!.......................................................................
!   (following comments added by rc_easter)
!
!   computes the following coagulation rate "coefficients"
!   for intramodal coagulation
!	dNX/dt  = - (1.e-20*quadn11)*NX1*NX1
!	dM2X/dt = - (1.e-20*quads11)*NX1*NX1
!
!   NX  is the mode's number concentration (particles/m3)
!   M2X is the mode's 2nd moment (surface/pi) concentration (m2/m3)
!
!   input arguments
!	lamda = mean free path (m)
!	kfm, knc = constants used in free-molecular and near-continuum
!	    calculations (see subr bink_coag_rates)
!	dg = DRY median diameter for number distribution (m)
!	xlnsig = ln of geometric standard deviation (dimensionless)
!	wetddryd = (WET median diameter)/(DRY median diameter)
!
!   output arguments
!	quads11, quadn11 = coagulation rate coefficients (see above)
!
!   *** note that the quads11, quadn11
!	are all scaled by 1.e+20 to avoid underflow, and they
!	must be multiplied by 1.e-20 when they are applied
!
!   *** wetddryd1 was added because MIRAGE treats N/S/V of
!	the DRY size distribution, so the quadrature should be done
!	over the DRY size distribution.  However, the coagulation kernel
!	must be computed using actual (WET) particle sizes
!
!.......................................................................
!
      implicit none

      integer i,j

      real(r4) lamda ! mean free path
      real(r4) kfm, knc
      real(r4) dg, xlnsig, wetddryd
      real(r4) quads11, quadn11
      real(r4) pi
      parameter( pi = 3.14159265358979)
      real(r4) two3rds
      parameter( two3rds = 2.0d0 /  3.0d0 )
      real(r4) sqrt2
      parameter(sqrt2 = 1.41421356237309 )
      real(r4) sum1sfm, sum2sfm, sum1nfm, sum2nfm
      real(r4) sum1snc, sum2snc, sum1nnc, sum2nnc
      real(r4) xi, wxi, xf, dp1p,dp1m,dp1psq,dp1msq, dp1pwet, dp1mwet
      real(r4) v1p,v1m, a2p,a2m,v2p,v2m
      real(r4) yi,wyi,yf,dp2p,dp2m,dp2psq,dp2msq, dp2pwet, dp2mwet
      real(r4) dspp,dsmp,dspm, dsmm
      real(r4) bppfm,bmpfm,bpmfm,bmmfm
      real(r4) bppnc,bmpnc,bpmnc,bmmnc
      real(r4) xx1, xx2
      real(r4) xbsfm, xbsnc, xbnfm, xbnnc
      real(r4) betafm, betanc

      real(r4)   a               ! approx cunningham corr. factor
      parameter( a = 1.246d0 )

      real(r4)   twoa
      parameter( twoa = 2.0d0 * a )

! *** has a fixed number of gauss-herimite abscissas (n)
      integer n    ! one-half the number of abscissas
      parameter ( n = 5 )
      real(r4), save :: ghxi(n) ! gauss-hermite abscissas
      real(r4), save :: ghwi(n) ! gauss-hermite weights

! ** values from table 25.10 (page 924) of abramowitz and stegun,
!    handbook of mathematical functions, national bureau of standards,
!  december 1965.

!    breaks in number to facilitate comparison with printed table

! *** tests show that 10 point is adquate.

      data ghxi/0.342901327223705,   &
                1.036610829789514,   &
                1.756683649299882,   &
                2.532731674232790,   &
                3.436159118837738/

      data ghwi/6.108626337353d-1,   &
                2.401386110823d-1,   &
                3.387439445548d-2,   &
                1.343645746781d-3,   &
                7.640432855233d-6/

! *** the following expressions are from binkowski & shanker
!     jour. geophys. research. vol. 100,no. d12, pp 26,191-26,209
!     december 20, 1995
! ***  for free molecular eq. a5
        betafm(xx1, xx2) = kfm *   &
             sqrt(1.0 / xx1**3  + 1.0 / xx2**3 ) * (xx1 + xx2)**2

! ***  for near continuum  eq. a6
        betanc(xx1, xx2) =  knc * (xx1 + xx2) *   &
                             ( 1.0d0 / xx1 + 1.0d0 / xx2  +   &
                           twoa * lamda * ( 1.0d0 / xx1 ** 2   &
                                          + 1.0d0 / xx2 **2 ) )


      sum1sfm = 0.0
      sum1snc = 0.0

      sum1nfm = 0.0
      sum1nnc = 0.0
      do 1 i=1,n

        sum2sfm = 0.0
        sum2snc = 0.0
        sum2nfm = 0.0
        sum2nnc = 0.0

        xi = ghxi(i)
        wxi = ghwi(i)
        xf = exp( sqrt2 * xi *xlnsig)
        dp1p = dg*xf
        dp1m = dg/xf
        dp1psq = dp1p*dp1p
        dp1msq = dp1m*dp1m
        v1p = dp1p*dp1psq
        v1m = dp1m*dp1msq

        dp1pwet = dp1p * wetddryd
        dp1mwet = dp1m * wetddryd

      do 11 j=1,n
        yi = ghxi(j)
        wyi = ghwi(j)
        yf = exp( sqrt2 * yi * xlnsig)
        dp2p = dg*yf
        dp2m = dg/yf
        dp2psq = dp2p*dp2p
        dp2msq = dp2m*dp2m
        a2p = dp2psq
        a2m = dp2msq
        v2p =  dp2p*dp2psq
        v2m =dp2m*dp2msq
        dspp = 0.5*(v1p+v2p)**two3rds - a2p
        dsmp = 0.5*(v1m+v2p)**two3rds - a2p
        dspm = 0.5*(v1p+v2m)**two3rds - a2m
        dsmm = 0.5*(v1m+v2m)**two3rds - a2m

        dp2pwet = dp2p * wetddryd
        dp2mwet = dp2m * wetddryd

!   scale by 1.0e+20 to avoid underflow
        bppfm = betafm(dp1pwet,dp2pwet) * 1.0e20
        bmpfm = betafm(dp1mwet,dp2pwet) * 1.0e20
        bpmfm = betafm(dp1pwet,dp2mwet) * 1.0e20
        bmmfm = betafm(dp1mwet,dp2mwet) * 1.0e20

        bppnc = betanc(dp1pwet,dp2pwet) * 1.0e20
        bmpnc = betanc(dp1mwet,dp2pwet) * 1.0e20
        bpmnc = betanc(dp1pwet,dp2mwet) * 1.0e20
        bmmnc = betanc(dp1mwet,dp2mwet) * 1.0e20

        sum2sfm = sum2sfm + wyi*(dspp * bppfm + dspm * bpmfm   &
                     +   dsmp * bmpfm + dsmm * bmmfm )

        sum2nfm = sum2nfm + wyi*(bppfm + bmpfm + bpmfm + bmmfm)

        sum2snc = sum2snc + wyi*(dspp * bppnc + dspm * bpmnc   &
                     +   dsmp * bmpnc + dsmm * bmmnc )
        sum2nnc = sum2nnc + wyi*(bppnc + bmpnc + bpmnc + bmmnc)

   11 continue
      sum1sfm = sum1sfm + wxi * sum2sfm
      sum1nfm = sum1nfm + wxi * sum2nfm

      sum1snc = sum1snc + wxi * sum2snc
      sum1nnc = sum1nnc + wxi * sum2nnc

    1 continue

      xbsfm   = -sum1sfm  / pi
      xbsnc   = -sum1snc  / pi

!     quads11 =  xbsfm * xbsnc / ( xbsfm + xbsnc )
      quads11 =  ( max(xbsfm,xbsnc) / ( xbsfm + xbsnc ) )   &
                 * min(xbsfm,xbsnc)

! *** quads11 is the intra-modal coagulation term for 2nd moment

      xbnfm   = 0.5 * sum1nfm  / pi
      xbnnc   = 0.5 * sum1nnc  / pi


!     quadn11 =  xbnfm * xbnnc / ( xbnfm + xbnnc )
      quadn11 =  ( max(xbnfm,xbnnc) / ( xbnfm + xbnnc ) )   &
                 * min(xbnfm,xbnnc)

! *** quadn11 is the intra-modal coagulation term for number


      return
      end subroutine intracoag_gh


! ---------------------------------------------------------------------
       subroutine intercoag_gh( lamda, kfm, knc, dg1, dg2,   &
           xlnsig1, xlnsig2,   &
           wetddryd1, wetddryd2,   &
           quads12, quads21, quadn12, quadv12 )

!  fsb this version is for intermodal coagulation for number,
!       second, and third moments
! fsb this version runs in real*4 arithmetic

! *** this version calculates the coagulation coefficients
!     using the harmonic mean approach for both fm and nc cases.
! *** does gauss-hermite quadrature for inter-modal
!      coagulation integrals for 2nd moment
!     for two lognormal distributions
!      defined by dg1,xlnsig1, dg2,xlnsig2
! *** dg and xlnsig are the geometric mean diameters (meters)
!      and the logarithms of the
!     geometric standard deviations (dimensionless)
!     whose meaning is defined below at the end of the routine
!     ghxi, ghwi are the gauss-hermite weights and n is one-half the
!     number of abscissas, since an even number of abscissas is used
!
!.......................................................................
!   (following comments added by rc_easter)
!
!   computes the following coagulation rate "coefficients"
!   for mode1-mode2 coagulation
!	dNX1/dt  = - (1.e-20*quadn12)*NX1*NX2
!	dNX2/dt  = 0.
!	dM2X1/dt = - (1.e-20*quads12)*NX1*NX2
!	dM2X2/dt = + (1.e-20*quads21)*NX1*NX2
!	dM3X1/dt = - (1.e-20*quads12)*NX1*NX2 == - dM3X2/dt
!
!   NX1,  NX2 are number concentrations (particles/m3)
!   M2X1, M2X2 are 2nd moment (surface/pi)  concentrations (m2/m3)
!   M3X1, M3X2 are 2rd moment (volume*6/pi) concentrations (m3/m3)
!
!   in the above, mode 1 is aitken mode, mode 2 is accumulation mode
!
!   input arguments
!	lamda = mean free path (m)
!	kfm, knc = constants used in free-molecular and near-continuum
!	    calculations (see subr bink_coag_rates)
!	dg1, dg2 = DRY median diameters for number distribution (m)
!	    for aitken and accumulation modes
!	xlnsig1, xlnsig2 = ln of geometric standard deviations (dimensionless)
!	wetddryd1, wetddryd2 = (WET median diameter)/(DRY median diameter)
!
!   output arguments
!	quads12, quads21, quadn12, quadv12 = coagulation rate
!	    coefficients (see above)
!
!   *** note that the quads12, quads21, quadn12, quadv12
!	are all scaled by 1.e+20 to avoid underflow, and they
!	must be multiplied by 1.e-20 when they are applied
!
!   *** wetddryd1/2 were added because MIRAGE treats N/S/V of
!	the DRY size distribution, so the quadrature should be done
!	over the DRY size distribution.  However, the coagulation kernel
!	must be computed using actual (WET) particle sizes
!
!.......................................................................

      implicit none

      integer i,j

      real(r4) lamda ! mean free path
      real(r4) kfm, knc
      real(r4) dg1, xlnsig1, dg2, xlnsig2
      real(r4) wetddryd1, wetddryd2
      real(r4) quads12, quads21, quadn12, quadv12
      real(r4) pi
      parameter( pi = 3.14159265358979)
      real(r4) two3rds
      parameter( two3rds = 2.0d0 /  3.0d0 )
      real(r4) sqrt2
      parameter(sqrt2 = 1.41421356237309 )
      real(r4) sum1s12fm, sum1s21fm, sum2s12fm, sum2s21fm
      real(r4) sum1nfm, sum2nfm
      real(r4) sum1vfm, sum2vfm
      real(r4) sum1s12nc, sum1s21nc, sum2s12nc, sum2s21nc
      real(r4) sum1nnc, sum2nnc
      real(r4) sum1vnc, sum2vnc
      real(r4) xi, wxi,xf, dp1p, dp1m, dp1psq, dp1msq, dp1pwet, dp1mwet
      real(r4) a1p, a1m, v1p, v1m
      real(r4) a2p, a2m, v2p, v2m
      real(r4) yi, wyi, yf, dp2p, dp2m, dp2psq, dp2msq, dp2pwet, dp2mwet
      real(r4) dspp, dsmp, dspm, dsmm
      real(r4) bppfm, bmpfm, bpmfm, bmmfm
      real(r4) bppnc, bmpnc, bpmnc, bmmnc
      real(r4) xx1, xx2
      real(r4) xbsfm, xbsnc, xbnfm, xbnnc, xbvfm, xbvnc
      real(r4) betafm, betanc

      real(r4)   a               ! approx cunningham corr. factor
      parameter( a = 1.246d0 )

      real(r4)   twoa
      parameter( twoa = 2.0d0 * a )

! *** has a fixed number of gauss-herimite abscissas (n)
      integer n    ! one-half the number of abscissas
      parameter ( n = 5 )
      real(r4), save :: ghxi(n) ! gauss-hermite abscissas
      real(r4), save :: ghwi(n) ! gauss-hermite weights

! ** values from table 25.10 (page 924) of abramowitz and stegun,
!    handbook of mathematical functions, national bureau of standards,
!  december 1965.

!    breaks in number to facilitate comparison with printed table

! *** tests show that 10 point is adquate.

      data ghxi/0.342901327223705,   &
                1.036610829789514,   &
                1.756683649299882,   &
                2.532731674232790,   &
                3.436159118837738/    

      data ghwi/6.108626337353d-1,   &
                2.401386110823d-1,   &
                3.387439445548d-2,   &
                1.343645746781d-3,   &
                7.640432855233d-6/


! *** the following expressions are from binkowski & shanker
!     jour. geophys. research. vol. 100,no. d12, pp 26,191-26,209
!     december 20, 1995

! ***  for free molecular eq. a5
        betafm(xx1, xx2) = kfm *   &
             sqrt(1.0 / xx1**3  + 1.0 / xx2**3 ) * (xx1 + xx2)**2

! ***  for near continuum  eq. a6
        betanc(xx1, xx2) =  knc * (xx1 + xx2) *   &
                             ( 1.0d0 / xx1 + 1.0d0 / xx2  +   &
                           twoa * lamda * ( 1.0d0 / xx1 ** 2   &
                                          + 1.0d0 / xx2 **2 ) )

      sum1s12fm = 0.0
      sum1s12nc = 0.0
      sum1s21fm = 0.0
      sum1s21nc = 0.0
      sum1vnc = 0.0
      sum1vfm = 0.0
      sum1nfm = 0.0
      sum1nnc = 0.0
      do 1 i=1,n

        sum2s12fm = 0.0
        sum2s12nc = 0.0
        sum2s21fm = 0.0
        sum2s21nc = 0.0
        sum2nfm = 0.0
        sum2nnc = 0.0
        sum2vnc = 0.0
        sum2vfm = 0.0
        xi = ghxi(i)
        wxi = ghwi(i)
        xf = exp( sqrt2 * xi *xlnsig1)
        dp1p = dg1*xf
        dp1m = dg1/xf
        dp1psq = dp1p*dp1p
        dp1msq = dp1m*dp1m
        a1p = dp1psq
        a1m = dp1msq
        v1p = dp1p*dp1psq
        v1m = dp1m*dp1msq

        dp1pwet = dp1p * wetddryd1
        dp1mwet = dp1m * wetddryd1

      do 11 j=1,n
        yi  = ghxi(j)
        wyi = ghwi(j)
        yf = exp( sqrt2 * yi * xlnsig2)
        dp2p = dg2*yf
        dp2m = dg2/yf
        dp2psq = dp2p*dp2p
        dp2msq = dp2m*dp2m
        a2p  = dp2psq
        a2m  = dp2msq
        v2p  =  dp2p*dp2psq
        v2m  = dp2m*dp2msq
        dspp = (v1p+v2p)**two3rds - a2p
        dsmp = (v1m+v2p)**two3rds - a2p
        dspm = (v1p+v2m)**two3rds - a2m
        dsmm = (v1m+v2m)**two3rds - a2m

        dp2pwet = dp2p * wetddryd2
        dp2mwet = dp2m * wetddryd2

!   scale by 1.0e+20 to avoid underflow
        bppfm = betafm(dp1pwet,dp2pwet) * 1.0e20
        bmpfm = betafm(dp1mwet,dp2pwet) * 1.0e20
        bpmfm = betafm(dp1pwet,dp2mwet) * 1.0e20
        bmmfm = betafm(dp1mwet,dp2mwet) * 1.0e20

        bppnc = betanc(dp1pwet,dp2pwet) * 1.0e20
        bmpnc = betanc(dp1mwet,dp2pwet) * 1.0e20
        bpmnc = betanc(dp1pwet,dp2mwet) * 1.0e20
        bmmnc = betanc(dp1mwet,dp2mwet) * 1.0e20


        sum2s12fm = sum2s12fm + wyi*(a1p * bppfm + a1p * bpmfm   &
                     +   a1m * bmpfm + a1m * bmmfm )

        sum2s21fm = sum2s21fm + wyi*(dspp * bppfm + dspm * bpmfm   &
                     +   dsmp * bmpfm + dsmm * bmmfm )


        sum2s12nc = sum2s12nc + wyi*(a1p * bppnc + a1p * bpmnc   &
                     +   a1m * bmpnc + a1m * bmmnc )

        sum2s21nc = sum2s21nc + wyi*(dspp * bppnc + dspm * bpmnc   &
                     +   dsmp * bmpnc + dsmm * bmmnc )

        sum2nfm = sum2nfm + wyi*(bppfm + bmpfm + bpmfm + bmmfm)

        sum2nnc = sum2nnc + wyi*(bppnc + bmpnc + bpmnc + bmmnc)

        sum2vfm = sum2vfm + wyi*(v1p*(bppfm + bpmfm) +   &
                                 v1m*(bmpfm + bmmfm) )

        sum2vnc = sum2vnc + wyi*(v1p*(bppnc + bpmnc) +   &
                                 v1m*(bmpnc + bmmnc) )

   11 continue

      sum1s12fm = sum1s12fm + wxi * sum2s12fm
      sum1s21fm = sum1s21fm + wxi * sum2s21fm
      sum1nfm   = sum1nfm + wxi * sum2nfm
      sum1vfm   = sum1vfm + wxi * sum2vfm

      sum1s12nc = sum1s12nc + wxi * sum2s12nc
      sum1s21nc = sum1s21nc + wxi * sum2s21nc
      sum1nnc   = sum1nnc + wxi * sum2nnc
      sum1vnc   = sum1vnc + wxi * sum2vnc

    1 continue

! *** second moment intermodal coagulation coefficients

! fsb note: the transfer of second moment is not symmetric.
!     see equations a3 & a4 of binkowski & shankar (1995)

! ***  to accumulation mode from aitken mode

      xbsfm   = sum1s21fm  / pi
      xbsnc   = sum1s21nc  / pi

!     quads21 =  xbsfm * xbsnc / ( xbsfm + xbsnc )
      quads21 =  ( max(xbsfm,xbsnc) / ( xbsfm + xbsnc ) )   &
                 * min(xbsfm,xbsnc)

! *** from aitken mode to accumulation mode

      xbsfm   = sum1s12fm  / pi
      xbsnc   = sum1s12nc  / pi

!     quads12 =  xbsfm * xbsnc / ( xbsfm + xbsnc )
      quads12 =  ( max(xbsfm,xbsnc) / ( xbsfm + xbsnc ) )   &
                 * min(xbsfm,xbsnc)



      xbnfm   = sum1nfm  / pi
      xbnnc   = sum1nnc  / pi

!     quadn12 =  xbnfm * xbnnc / ( xbnfm + xbnnc )
      quadn12 =  ( max(xbnfm,xbnnc) / ( xbnfm + xbnnc ) )   &
                 * min(xbnfm,xbnnc)

! *** quadn12 is the intermodal coagulation coefficient for number


       xbvfm = sum1vfm / pi
       xbvnc = sum1vnc / pi

!     quadv12 = xbvfm * xbvnc / ( xbvfm + xbvnc )
      quadv12 =  ( max(xbvfm,xbvnc) / ( xbvfm + xbvnc ) )   &
                 * min(xbvfm,xbvnc)

! *** quadv12 is the intermodal coagulation coefficient for 3rd moment


      return
      end subroutine intercoag_gh


!----------------------------------------------------------------------
!----------------------------------------------------------------------
      subroutine getcoags_wrapper_f(              &
          airtemp, airprs,                        &
          dgatk, dgacc,                           &
          sgatk, sgacc,                           &
          xxlsgat, xxlsgac,                       &
          pdensat, pdensac,                       &
          betaij0, betaij2i, betaij2j, betaij3,   &
          betaii0, betaii2, betajj0, betajj2      )
!
! interface to subr. getcoags
!
! interface code adapted from subr. aeroproc of cmaq v4.6,
!     with some of the parameter values from module aero_info_ae4
!
      implicit none

! *** arguments

      real(r8), intent(in) :: airtemp  ! air temperature [ k ]
      real(r8), intent(in) :: airprs   ! air pressure in [ pa ]

      real(r8), intent(in) :: dgatk    ! aitken mode geometric mean diameter [m]
      real(r8), intent(in) :: dgacc    ! accumulation mode geometric mean diam [m]

      real(r8), intent(in) :: sgatk    ! aitken mode geometric standard deviation
      real(r8), intent(in) :: sgacc    ! accumulation mode geometric standard deviation

      real(r8), intent(in) :: xxlsgat  ! natural log of geometric standard
      real(r8), intent(in) :: xxlsgac  !  deviations

      real(r8), intent(in) :: pdensat  ! aitken mode particle density [ kg / m**3 ]
      real(r8), intent(in) :: pdensac  ! accumulation mode density [ kg / m**3 ]

      real(r8), intent(out) :: betaij0, betaij2i, betaij2j, betaij3,   &
                               betaii0, betaii2,  betajj0,  betajj2


! *** local parameters
      real(r8), parameter :: p0 = 101325.0         !  standard surface pressure [ pa ]
      real(r8), parameter :: t0 = 288.15           !  standard surface temperature [ k ]
      real(r8), parameter :: avo = 6.0221367e23    ! avogadro's constant [ 1/mol ]
      real(r8), parameter :: rgasuniv = 8.314510   ! universal gas constant [ j/mol-k ]
      real(r8), parameter :: boltz = rgasuniv/avo  ! boltzmann's constant [ j / k]
      real(r8), parameter :: two3 = 2.0/3.0

! *** local variables
      real(r8) amu            ! atmospheric dynamic viscosity [ kg/m s ]
      real(r8) sqrt_temp      ! square root of ambient temperature
      real(r8) lamda          ! mean free path [ m ]

! *** intramodal coagulation rates [ m**3/s ] ( 0th & 2nd moments )
      real(r8)    batat( 2 )  ! aitken mode
      real(r8)    bacac( 2 )  ! accumulation mode
! *** intermodal coagulation rates [ m**3/s ] ( 0th & 2nd moments )
      real(r8)    batac( 2 )  ! aitken to accumulation
      real(r8)    bacat( 2 )  ! accumulation from aitken
! *** intermodal coagulation rate [ m**3/s ] ( 3rd moment )
      real(r8)    c3ij        ! aitken to accumulation
! *** 3rd moment intermodal transfer rate by coagulation
      real(r8)    c30atac     ! aitken to accumulation

! *** near continnuum regime (independent of mode)
      real(r8)    knc         ! knc = two3 * boltz *  airtemp / amu
! *** free molecular regime (depends upon modal density)
      real(r8)    kfmat       ! kfmat = sqrt(3.0*boltz*airtemp/pdensat)
      real(r8)    kfmac       ! kfmac = sqrt(3.0*boltz*airtemp/pdensac)
      real(r8)    kfmatac     ! kfmatac = sqrt( 6.0 * boltz * airtemp /
                              !                ( pdensat + pdensac ) )

      real(r8)    dumacc2, dumatk2, dumatk3



      sqrt_temp = sqrt( airtemp)

! *** calculate mean free path [ m ]:
!     6.6328e-8 is the sea level value given in table i.2.8
!     on page 10 of u.s. standard atmosphere 1962
      lamda = 6.6328e-8 * p0 * airtemp  / ( t0 * airprs )

! *** calculate dynamic viscosity [ kg m**-1 s**-1 ]:
!     u.s. standard atmosphere 1962 page 14 expression
!     for dynamic viscosity is:
!     dynamic viscosity =  beta * t * sqrt(t) / ( t + s)
!     where beta = 1.458e-6 [ kg sec^-1 k**-0.5 ], s = 110.4 [ k ].
      amu = 1.458e-6 * airtemp * sqrt_temp / ( airtemp + 110.4 )

! *** coagulation
!     calculate coagulation coefficients using a method dictated by
!     the value of fastcoag_flag.  if true, the computationally-
!     efficient getcoags routine is used.  if false, the more intensive
!     gauss-hermite numerical quadrature method is used.  see section
!     2.1 of bhave et al. (2004) for further discussion.

! *** calculate term used in equation a6 of binkowski & shankar (1995)
      knc      = two3 * boltz *  airtemp / amu
! *** calculate terms used in equation a5 of binkowski & shankar (1995)
      kfmat    = sqrt( 3.0 * boltz * airtemp / pdensat )
      kfmac    = sqrt( 3.0 * boltz * airtemp / pdensac )
      kfmatac  = sqrt( 6.0 * boltz * airtemp / ( pdensat + pdensac ) )

! *** transfer of number to accumulation mode from aitken mode is zero
      bacat(1) = 0.0

! *** calculate intermodal and intramodal coagulation coefficients
!     for zeroth and second moments, and intermodal coagulation
!     coefficient for third moment
        call getcoags( lamda, kfmatac, kfmat, kfmac, knc,   &
                       dgatk,   dgacc,   sgatk,   sgacc,     &
                       xxlsgat,  xxlsgac,     &
                       batat(2), batat(1), bacac(2), bacac(1),   &
                       batac(2), bacat(2), batac(1), c3ij )

! convert from the "cmaq" coag rate parameters 
! to the "mirage2" parameters
        dumacc2 = ( (dgacc**2) * exp( 2.0*xxlsgac*xxlsgac ) )
        dumatk2 = ( (dgatk**2) * exp( 2.0*xxlsgat*xxlsgat ) )
        dumatk3 = ( (dgatk**3) * exp( 4.5*xxlsgat*xxlsgat ) )

        betaii0  = max( 0.0_r8, batat(1) )
        betajj0  = max( 0.0_r8, bacac(1) )
        betaij0  = max( 0.0_r8, batac(1) )
        betaij3  = max( 0.0_r8, c3ij / dumatk3 )

        betajj2  = max( 0.0_r8, bacac(2) / dumacc2 )
        betaii2  = max( 0.0_r8, batat(2) / dumatk2 )
        betaij2i = max( 0.0_r8, batac(2) / dumatk2 )
        betaij2j = max( 0.0_r8, bacat(2) / dumatk2 )


      return
      end subroutine getcoags_wrapper_f



! //////////////////////////////////////////////////////////////////
!  subroutine getcoags calculates the coagulation rates using a new
!     approximate algorithm for the 2nd moment.  the 0th and 3rd moments
!     are done by analytic expressions from whitby et al. (1991).  the
!     correction factors are also similar to those from whitby et al.
!     (1991), but are derived from the gauss-hermite numerical
!     quadratures used by binkowski and roselle (2003).
!
!     called from aerostep as:
!     call getcoags( lamda, kfmatac, kfmat, kfmac, knc,
!                    dgat,dgac, sgatk, sgacc, xxlsgat,xxlsgac,
!                    batat(2), batat(1), bacac(2), bacac(1),
!                    batac(2), bacat(2), batac(1), c3ij )
!     where all input and outputs are real*8
!
!  revision history:
!   fsb 08/25/03 coded by dr. francis s. binkowksi
!
!   fsb 08/25/04 added in-line documentation
!
!   rce 04/15/2007 
!       code taken from cmaq v4.6 code; converted to f90;
!	added "intent" to subr arguments;
!       renamed "r4" & "r8" variables to "rx4" & "rx8";
!       changed "real*N" declarations to "real(rN)" (N = 4 or 8)
!
!  references:
!   1. whitby, e. r., p. h. mcmurry, u. shankar, and f. s. binkowski,
!   modal aerosol dynamics modeling, rep. 600/3-91/020, atmospheric
!   research and exposure assessment laboratory,
!   u.s. environmental protection agency, research triangle park, n.c.,
!   (ntis pb91-161729/as), 1991
!
!   2. binkowski, f.s. an u. shankar, the regional particulate matter
!   model 1. model decsription and preliminary results, journal of
!   geophysical research, 100, d12, pp 26,191-26,209,
!   december 20, 1995.
!
!   3. binkowski, f.s. and s.j. roselle, models-3 community
!      multiscale air quality (cmaq) model aerosol component 1:
!      model description.  j. geophys. res., vol 108, no d6, 4183
!      doi:10.1029/2001jd001409, 2003.


      subroutine getcoags( lamda, kfmatac, kfmat, kfmac, knc,   &
                           dgatk, dgacc, sgatk, sgacc, xxlsgat,xxlsgac,   &
                           qs11, qn11, qs22, qn22,   &
                           qs12, qs21, qn12, qv12 )

      implicit none

      real(r8), intent(in) ::  lamda     ! mean free path [ m ]

! *** coefficients for free molecular regime
      real(r8), intent(in) ::  kfmat     ! aitken mode
      real(r8), intent(in) ::  kfmac     ! accumulation mode
      real(r8), intent(in) ::  kfmatac   ! aitken to accumulation mode

      real(r8), intent(in) ::  knc   ! coefficient for near continnuum regime

! *** modal geometric mean diameters: [ m ]
      real(r8), intent(in) :: dgatk          ! aitken mode
      real(r8), intent(in) :: dgacc          ! accumulation mode

! *** modal geometric standard deviation
      real(r8), intent(in) :: sgatk          ! atken mode
      real(r8), intent(in) :: sgacc          ! accumulation mode

! *** natural log of modal geometric standard deviation
      real(r8), intent(in) :: xxlsgat         ! aitken mode
      real(r8), intent(in) :: xxlsgac         ! accumulation mode

! *** coagulation coefficients
      real(r8), intent(out) :: qs11, qn11, qs22, qn22,   &
                               qs12, qs21, qn12, qv12

      integer ibeta, n1, n2a, n2n ! indices for correction factors

      real(r8)  i1fm_at
      real(r8)  i1nc_at
      real(r8)  i1_at

      real(r8)  i1fm_ac
      real(r8)  i1nc_ac
      real(r8)  i1_ac

      real(r8)  i1fm
      real(r8)  i1nc
      real(r8)  i1

      real(r8) constii

      real(r8)    kngat, kngac
      real(r8)    one, two, half
       parameter( one = 1.0d0, two = 2.0d0, half = 0.5d0 )
      real(r8)    a
!       parameter( a = 2.492d0)
      parameter( a = 1.246d0)
      real      two3rds
       parameter( two3rds = 2.d0 / 3.d0)

      real(r8)   sqrttwo  !  sqrt(two)
      real(r8)   dlgsqt2  !  1/ln( sqrt( 2 ) )


      real(r8)    esat01         ! aitken mode exp( log^2( sigmag )/8 )
      real(r8)    esac01         ! accumulation mode exp( log^2( sigmag )/8 )

      real(r8)    esat04
      real(r8)    esac04

      real(r8)    esat05
      real(r8)    esac05

      real(r8)    esat08
      real(r8)    esac08

      real(r8)    esat09
      real(r8)    esac09

      real(r8)    esat16
      real(r8)    esac16

      real(r8)    esat20
      real(r8)    esac20

      real(r8)    esat24
      real(r8)    esac24

      real(r8)    esat25
      real(r8)    esac25

      real(r8)    esat36
      real(r8)    esac36

      real(r8)    esat49

      real(r8)    esat64
      real(r8)    esac64

      real(r8)    esat100

      real(r8) dgat2, dgac2, dgat3, dgac3
      real(r8) sqdgat, sqdgac
      real(r8) sqdgat5, sqdgac5
      real(r8) sqdgat7
      real(r8) r, r2, r3, rx4, r5, r6, rx8
      real(r8) ri1, ri2, ri3, ri4
      real(r8) rat
      real(r8) coagfm0, coagnc0
      real(r8) coagfm3, coagnc3
      real(r8) coagfm_at, coagfm_ac
      real(r8) coagnc_at, coagnc_ac
      real(r8) coagatat0
      real(r8) coagacac0
      real(r8) coagatat2
      real(r8) coagacac2
      real(r8) coagatac0, coagatac3
      real(r8) coagatac2
      real(r8) coagacat2
      real(r8) xm2at, xm3at, xm2ac, xm3ac

! *** correction factors for coagulation rates
      real(r4), save :: bm0( 10 )          ! m0 intramodal fm - rpm values
      real(r4), save :: bm0ij( 10, 10, 10 ) ! m0 intermodal fm
      real(r4), save :: bm3i( 10, 10, 10 ) ! m3 intermodal fm- rpm values
      real(r4), save :: bm2ii(10) ! m2 intramodal fm
      real(r4), save :: bm2iitt(10) ! m2 intramodal total
      real(r4), save :: bm2ij(10,10,10) ! m2 intermodal fm i to j
      real(r4), save :: bm2ji(10,10,10) ! m2 total intermodal  j from i

! *** populate the arrays for the correction factors.

! rpm 0th moment correction factors for unimodal fm coagulation  rates
      data      bm0  /   &
            0.707106785165097, 0.726148960080488, 0.766430744110958,   &
            0.814106389441342, 0.861679526483207, 0.903600509090092,   &
            0.936578814219156, 0.960098926735545, 0.975646823342881,   &
            0.985397173215326   /


! fsb new fm correction factors for m0 intermodal coagulation

      data (bm0ij (  1,  1,ibeta), ibeta = 1,10) /   &
        0.628539,  0.639610,  0.664514,  0.696278,  0.731558,   &
        0.768211,  0.804480,  0.838830,  0.870024,  0.897248/
      data (bm0ij (  1,  2,ibeta), ibeta = 1,10) /   &
        0.639178,  0.649966,  0.674432,  0.705794,  0.740642,   &
        0.776751,  0.812323,  0.845827,  0.876076,  0.902324/
      data (bm0ij (  1,  3,ibeta), ibeta = 1,10) /   &
        0.663109,  0.673464,  0.697147,  0.727637,  0.761425,   &
        0.796155,  0.829978,  0.861419,  0.889424,  0.913417/
      data (bm0ij (  1,  4,ibeta), ibeta = 1,10) /   &
        0.693693,  0.703654,  0.726478,  0.755786,  0.787980,   &
        0.820626,  0.851898,  0.880459,  0.905465,  0.926552/
      data (bm0ij (  1,  5,ibeta), ibeta = 1,10) /   &
        0.727803,  0.737349,  0.759140,  0.786870,  0.816901,   &
        0.846813,  0.874906,  0.900060,  0.921679,  0.939614/
      data (bm0ij (  1,  6,ibeta), ibeta = 1,10) /   &
        0.763461,  0.772483,  0.792930,  0.818599,  0.845905,   &
        0.872550,  0.897051,  0.918552,  0.936701,  0.951528/
      data (bm0ij (  1,  7,ibeta), ibeta = 1,10) /   &
        0.799021,  0.807365,  0.826094,  0.849230,  0.873358,   &
        0.896406,  0.917161,  0.935031,  0.949868,  0.961828/
      data (bm0ij (  1,  8,ibeta), ibeta = 1,10) /   &
        0.833004,  0.840514,  0.857192,  0.877446,  0.898147,   &
        0.917518,  0.934627,  0.949106,  0.960958,  0.970403/
      data (bm0ij (  1,  9,ibeta), ibeta = 1,10) /   &
        0.864172,  0.870734,  0.885153,  0.902373,  0.919640,   &
        0.935494,  0.949257,  0.960733,  0.970016,  0.977346/
      data (bm0ij (  1, 10,ibeta), ibeta = 1,10) /   &
        0.891658,  0.897227,  0.909343,  0.923588,  0.937629,   &
        0.950307,  0.961151,  0.970082,  0.977236,  0.982844/
      data (bm0ij (  2,  1,ibeta), ibeta = 1,10) /   &
        0.658724,  0.670587,  0.697539,  0.731890,  0.769467,   &
        0.807391,  0.843410,  0.875847,  0.903700,  0.926645/
      data (bm0ij (  2,  2,ibeta), ibeta = 1,10) /   &
        0.667070,  0.678820,  0.705538,  0.739591,  0.776758,   &
        0.814118,  0.849415,  0.881020,  0.908006,  0.930121/
      data (bm0ij (  2,  3,ibeta), ibeta = 1,10) /   &
        0.686356,  0.697839,  0.723997,  0.757285,  0.793389,   &
        0.829313,  0.862835,  0.892459,  0.917432,  0.937663/
      data (bm0ij (  2,  4,ibeta), ibeta = 1,10) /   &
        0.711425,  0.722572,  0.747941,  0.780055,  0.814518,   &
        0.848315,  0.879335,  0.906290,  0.928658,  0.946526/
      data (bm0ij (  2,  5,ibeta), ibeta = 1,10) /   &
        0.739575,  0.750307,  0.774633,  0.805138,  0.837408,   &
        0.868504,  0.896517,  0.920421,  0.939932,  0.955299/
      data (bm0ij (  2,  6,ibeta), ibeta = 1,10) /   &
        0.769143,  0.779346,  0.802314,  0.830752,  0.860333,   &
        0.888300,  0.913014,  0.933727,  0.950370,  0.963306/
      data (bm0ij (  2,  7,ibeta), ibeta = 1,10) /   &
        0.798900,  0.808431,  0.829700,  0.855653,  0.882163,   &
        0.906749,  0.928075,  0.945654,  0.959579,  0.970280/
      data (bm0ij (  2,  8,ibeta), ibeta = 1,10) /   &
        0.827826,  0.836542,  0.855808,  0.878954,  0.902174,   &
        0.923316,  0.941345,  0.955989,  0.967450,  0.976174/
      data (bm0ij (  2,  9,ibeta), ibeta = 1,10) /   &
        0.855068,  0.862856,  0.879900,  0.900068,  0.919956,   &
        0.937764,  0.952725,  0.964726,  0.974027,  0.981053/
      data (bm0ij (  2, 10,ibeta), ibeta = 1,10) /   &
        0.879961,  0.886755,  0.901484,  0.918665,  0.935346,   &
        0.950065,  0.962277,  0.971974,  0.979432,  0.985033/
      data (bm0ij (  3,  1,ibeta), ibeta = 1,10) /   &
        0.724166,  0.735474,  0.761359,  0.794045,  0.828702,   &
        0.862061,  0.891995,  0.917385,  0.937959,  0.954036/
      data (bm0ij (  3,  2,ibeta), ibeta = 1,10) /   &
        0.730416,  0.741780,  0.767647,  0.800116,  0.834344,   &
        0.867093,  0.896302,  0.920934,  0.940790,  0.956237/
      data (bm0ij (  3,  3,ibeta), ibeta = 1,10) /   &
        0.745327,  0.756664,  0.782255,  0.814026,  0.847107,   &
        0.878339,  0.905820,  0.928699,  0.946931,  0.960977/
      data (bm0ij (  3,  4,ibeta), ibeta = 1,10) /   &
        0.765195,  0.776312,  0.801216,  0.831758,  0.863079,   &
        0.892159,  0.917319,  0.937939,  0.954145,  0.966486/
      data (bm0ij (  3,  5,ibeta), ibeta = 1,10) /   &
        0.787632,  0.798347,  0.822165,  0.850985,  0.880049,   &
        0.906544,  0.929062,  0.947218,  0.961288,  0.971878/
      data (bm0ij (  3,  6,ibeta), ibeta = 1,10) /   &
        0.811024,  0.821179,  0.843557,  0.870247,  0.896694,   &
        0.920365,  0.940131,  0.955821,  0.967820,  0.976753/
      data (bm0ij (  3,  7,ibeta), ibeta = 1,10) /   &
        0.834254,  0.843709,  0.864356,  0.888619,  0.912245,   &
        0.933019,  0.950084,  0.963438,  0.973530,  0.980973/
      data (bm0ij (  3,  8,ibeta), ibeta = 1,10) /   &
        0.856531,  0.865176,  0.883881,  0.905544,  0.926290,   &
        0.944236,  0.958762,  0.969988,  0.978386,  0.984530/
      data (bm0ij (  3,  9,ibeta), ibeta = 1,10) /   &
        0.877307,  0.885070,  0.901716,  0.920729,  0.938663,   &
        0.953951,  0.966169,  0.975512,  0.982442,  0.987477/
      data (bm0ij (  3, 10,ibeta), ibeta = 1,10) /   &
        0.896234,  0.903082,  0.917645,  0.934069,  0.949354,   &
        0.962222,  0.972396,  0.980107,  0.985788,  0.989894/
      data (bm0ij (  4,  1,ibeta), ibeta = 1,10) /   &
        0.799294,  0.809144,  0.831293,  0.858395,  0.885897,   &
        0.911031,  0.932406,  0.949642,  0.963001,  0.973062/
      data (bm0ij (  4,  2,ibeta), ibeta = 1,10) /   &
        0.804239,  0.814102,  0.836169,  0.862984,  0.890003,   &
        0.914535,  0.935274,  0.951910,  0.964748,  0.974381/
      data (bm0ij (  4,  3,ibeta), ibeta = 1,10) /   &
        0.815910,  0.825708,  0.847403,  0.873389,  0.899185,   &
        0.922275,  0.941543,  0.956826,  0.968507,  0.977204/
      data (bm0ij (  4,  4,ibeta), ibeta = 1,10) /   &
        0.831348,  0.840892,  0.861793,  0.886428,  0.910463,   &
        0.931614,  0.948993,  0.962593,  0.972872,  0.980456/
      data (bm0ij (  4,  5,ibeta), ibeta = 1,10) /   &
        0.848597,  0.857693,  0.877402,  0.900265,  0.922180,   &
        0.941134,  0.956464,  0.968298,  0.977143,  0.983611/
      data (bm0ij (  4,  6,ibeta), ibeta = 1,10) /   &
        0.866271,  0.874764,  0.892984,  0.913796,  0.933407,   &
        0.950088,  0.963380,  0.973512,  0.981006,  0.986440/
      data (bm0ij (  4,  7,ibeta), ibeta = 1,10) /   &
        0.883430,  0.891216,  0.907762,  0.926388,  0.943660,   &
        0.958127,  0.969499,  0.978070,  0.984351,  0.988872/
      data (bm0ij (  4,  8,ibeta), ibeta = 1,10) /   &
        0.899483,  0.906505,  0.921294,  0.937719,  0.952729,   &
        0.965131,  0.974762,  0.981950,  0.987175,  0.990912/
      data (bm0ij (  4,  9,ibeta), ibeta = 1,10) /   &
        0.914096,  0.920337,  0.933373,  0.947677,  0.960579,   &
        0.971111,  0.979206,  0.985196,  0.989520,  0.992597/
      data (bm0ij (  4, 10,ibeta), ibeta = 1,10) /   &
        0.927122,  0.932597,  0.943952,  0.956277,  0.967268,   &
        0.976147,  0.982912,  0.987882,  0.991450,  0.993976/
      data (bm0ij (  5,  1,ibeta), ibeta = 1,10) /   &
        0.865049,  0.872851,  0.889900,  0.909907,  0.929290,   &
        0.946205,  0.959991,  0.970706,  0.978764,  0.984692/
      data (bm0ij (  5,  2,ibeta), ibeta = 1,10) /   &
        0.868989,  0.876713,  0.893538,  0.913173,  0.932080,   &
        0.948484,  0.961785,  0.972080,  0.979796,  0.985457/
      data (bm0ij (  5,  3,ibeta), ibeta = 1,10) /   &
        0.878010,  0.885524,  0.901756,  0.920464,  0.938235,   &
        0.953461,  0.965672,  0.975037,  0.982005,  0.987085/
      data (bm0ij (  5,  4,ibeta), ibeta = 1,10) /   &
        0.889534,  0.896698,  0.912012,  0.929395,  0.945647,   &
        0.959366,  0.970227,  0.978469,  0.984547,  0.988950/
      data (bm0ij (  5,  5,ibeta), ibeta = 1,10) /   &
        0.902033,  0.908713,  0.922848,  0.938648,  0.953186,   &
        0.965278,  0.974729,  0.981824,  0.987013,  0.990746/
      data (bm0ij (  5,  6,ibeta), ibeta = 1,10) /   &
        0.914496,  0.920599,  0.933389,  0.947485,  0.960262,   &
        0.970743,  0.978839,  0.984858,  0.989225,  0.992348/
      data (bm0ij (  5,  7,ibeta), ibeta = 1,10) /   &
        0.926281,  0.931761,  0.943142,  0.955526,  0.966600,   &
        0.975573,  0.982431,  0.987485,  0.991128,  0.993718/
      data (bm0ij (  5,  8,ibeta), ibeta = 1,10) /   &
        0.937029,  0.941877,  0.951868,  0.962615,  0.972112,   &
        0.979723,  0.985488,  0.989705,  0.992725,  0.994863/
      data (bm0ij (  5,  9,ibeta), ibeta = 1,10) /   &
        0.946580,  0.950819,  0.959494,  0.968732,  0.976811,   &
        0.983226,  0.988047,  0.991550,  0.994047,  0.995806/
      data (bm0ij (  5, 10,ibeta), ibeta = 1,10) /   &
        0.954909,  0.958581,  0.966049,  0.973933,  0.980766,   &
        0.986149,  0.990166,  0.993070,  0.995130,  0.996577/
      data (bm0ij (  6,  1,ibeta), ibeta = 1,10) /   &
        0.914182,  0.919824,  0.931832,  0.945387,  0.957999,   &
        0.968606,  0.976982,  0.983331,  0.988013,  0.991407/
      data (bm0ij (  6,  2,ibeta), ibeta = 1,10) /   &
        0.917139,  0.922665,  0.934395,  0.947580,  0.959792,   &
        0.970017,  0.978062,  0.984138,  0.988609,  0.991843/
      data (bm0ij (  6,  3,ibeta), ibeta = 1,10) /   &
        0.923742,  0.928990,  0.940064,  0.952396,  0.963699,   &
        0.973070,  0.980381,  0.985866,  0.989878,  0.992768/
      data (bm0ij (  6,  4,ibeta), ibeta = 1,10) /   &
        0.931870,  0.936743,  0.946941,  0.958162,  0.968318,   &
        0.976640,  0.983069,  0.987853,  0.991330,  0.993822/
      data (bm0ij (  6,  5,ibeta), ibeta = 1,10) /   &
        0.940376,  0.944807,  0.954004,  0.963999,  0.972928,   &
        0.980162,  0.985695,  0.989779,  0.992729,  0.994833/
      data (bm0ij (  6,  6,ibeta), ibeta = 1,10) /   &
        0.948597,  0.952555,  0.960703,  0.969454,  0.977181,   &
        0.983373,  0.988067,  0.991507,  0.993977,  0.995730/
      data (bm0ij (  6,  7,ibeta), ibeta = 1,10) /   &
        0.956167,  0.959648,  0.966763,  0.974326,  0.980933,   &
        0.986177,  0.990121,  0.992993,  0.995045,  0.996495/
      data (bm0ij (  6,  8,ibeta), ibeta = 1,10) /   &
        0.962913,  0.965937,  0.972080,  0.978552,  0.984153,   &
        0.988563,  0.991857,  0.994242,  0.995938,  0.997133/
      data (bm0ij (  6,  9,ibeta), ibeta = 1,10) /   &
        0.968787,  0.971391,  0.976651,  0.982148,  0.986869,   &
        0.990560,  0.993301,  0.995275,  0.996675,  0.997657/
      data (bm0ij (  6, 10,ibeta), ibeta = 1,10) /   &
        0.973822,  0.976047,  0.980523,  0.985170,  0.989134,   &
        0.992215,  0.994491,  0.996124,  0.997277,  0.998085/
      data (bm0ij (  7,  1,ibeta), ibeta = 1,10) /   &
        0.947410,  0.951207,  0.959119,  0.967781,  0.975592,   &
        0.981981,  0.986915,  0.990590,  0.993266,  0.995187/
      data (bm0ij (  7,  2,ibeta), ibeta = 1,10) /   &
        0.949477,  0.953161,  0.960824,  0.969187,  0.976702,   &
        0.982831,  0.987550,  0.991057,  0.993606,  0.995434/
      data (bm0ij (  7,  3,ibeta), ibeta = 1,10) /   &
        0.954008,  0.957438,  0.964537,  0.972232,  0.979095,   &
        0.984653,  0.988907,  0.992053,  0.994330,  0.995958/
      data (bm0ij (  7,  4,ibeta), ibeta = 1,10) /   &
        0.959431,  0.962539,  0.968935,  0.975808,  0.981882,   &
        0.986759,  0.990466,  0.993190,  0.995153,  0.996552/
      data (bm0ij (  7,  5,ibeta), ibeta = 1,10) /   &
        0.964932,  0.967693,  0.973342,  0.979355,  0.984620,   &
        0.988812,  0.991974,  0.994285,  0.995943,  0.997119/
      data (bm0ij (  7,  6,ibeta), ibeta = 1,10) /   &
        0.970101,  0.972517,  0.977428,  0.982612,  0.987110,   &
        0.990663,  0.993326,  0.995261,  0.996644,  0.997621/
      data (bm0ij (  7,  7,ibeta), ibeta = 1,10) /   &
        0.974746,  0.976834,  0.981055,  0.985475,  0.989280,   &
        0.992265,  0.994488,  0.996097,  0.997241,  0.998048/
      data (bm0ij (  7,  8,ibeta), ibeta = 1,10) /   &
        0.978804,  0.980591,  0.984187,  0.987927,  0.991124,   &
        0.993617,  0.995464,  0.996795,  0.997739,  0.998403/
      data (bm0ij (  7,  9,ibeta), ibeta = 1,10) /   &
        0.982280,  0.983799,  0.986844,  0.989991,  0.992667,   &
        0.994742,  0.996273,  0.997372,  0.998149,  0.998695/
      data (bm0ij (  7, 10,ibeta), ibeta = 1,10) /   &
        0.985218,  0.986503,  0.989071,  0.991711,  0.993945,   &
        0.995669,  0.996937,  0.997844,  0.998484,  0.998932/
      data (bm0ij (  8,  1,ibeta), ibeta = 1,10) /   &
        0.968507,  0.970935,  0.975916,  0.981248,  0.985947,   &
        0.989716,  0.992580,  0.994689,  0.996210,  0.997297/
      data (bm0ij (  8,  2,ibeta), ibeta = 1,10) /   &
        0.969870,  0.972210,  0.977002,  0.982119,  0.986619,   &
        0.990219,  0.992951,  0.994958,  0.996405,  0.997437/
      data (bm0ij (  8,  3,ibeta), ibeta = 1,10) /   &
        0.972820,  0.974963,  0.979339,  0.983988,  0.988054,   &
        0.991292,  0.993738,  0.995529,  0.996817,  0.997734/
      data (bm0ij (  8,  4,ibeta), ibeta = 1,10) /   &
        0.976280,  0.978186,  0.982060,  0.986151,  0.989706,   &
        0.992520,  0.994636,  0.996179,  0.997284,  0.998069/
      data (bm0ij (  8,  5,ibeta), ibeta = 1,10) /   &
        0.979711,  0.981372,  0.984735,  0.988263,  0.991309,   &
        0.993706,  0.995499,  0.996801,  0.997730,  0.998389/
      data (bm0ij (  8,  6,ibeta), ibeta = 1,10) /   &
        0.982863,  0.984292,  0.987172,  0.990174,  0.992750,   &
        0.994766,  0.996266,  0.997352,  0.998125,  0.998670/
      data (bm0ij (  8,  7,ibeta), ibeta = 1,10) /   &
        0.985642,  0.986858,  0.989301,  0.991834,  0.993994,   &
        0.995676,  0.996923,  0.997822,  0.998460,  0.998910/
      data (bm0ij (  8,  8,ibeta), ibeta = 1,10) /   &
        0.988029,  0.989058,  0.991116,  0.993240,  0.995043,   &
        0.996440,  0.997472,  0.998214,  0.998739,  0.999108/
      data (bm0ij (  8,  9,ibeta), ibeta = 1,10) /   &
        0.990046,  0.990912,  0.992640,  0.994415,  0.995914,   &
        0.997073,  0.997925,  0.998536,  0.998968,  0.999271/
      data (bm0ij (  8, 10,ibeta), ibeta = 1,10) /   &
        0.991732,  0.992459,  0.993906,  0.995386,  0.996633,   &
        0.997592,  0.998296,  0.998799,  0.999154,  0.999403/
      data (bm0ij (  9,  1,ibeta), ibeta = 1,10) /   &
        0.981392,  0.982893,  0.985938,  0.989146,  0.991928,   &
        0.994129,  0.995783,  0.996991,  0.997857,  0.998473/
      data (bm0ij (  9,  2,ibeta), ibeta = 1,10) /   &
        0.982254,  0.983693,  0.986608,  0.989673,  0.992328,   &
        0.994424,  0.995998,  0.997146,  0.997969,  0.998553/
      data (bm0ij (  9,  3,ibeta), ibeta = 1,10) /   &
        0.984104,  0.985407,  0.988040,  0.990798,  0.993178,   &
        0.995052,  0.996454,  0.997474,  0.998204,  0.998722/
      data (bm0ij (  9,  4,ibeta), ibeta = 1,10) /   &
        0.986243,  0.987386,  0.989687,  0.992087,  0.994149,   &
        0.995765,  0.996971,  0.997846,  0.998470,  0.998913/
      data (bm0ij (  9,  5,ibeta), ibeta = 1,10) /   &
        0.988332,  0.989313,  0.991284,  0.993332,  0.995082,   &
        0.996449,  0.997465,  0.998200,  0.998723,  0.999093/
      data (bm0ij (  9,  6,ibeta), ibeta = 1,10) /   &
        0.990220,  0.991053,  0.992721,  0.994445,  0.995914,   &
        0.997056,  0.997902,  0.998513,  0.998947,  0.999253/
      data (bm0ij (  9,  7,ibeta), ibeta = 1,10) /   &
        0.991859,  0.992561,  0.993961,  0.995403,  0.996626,   &
        0.997574,  0.998274,  0.998778,  0.999136,  0.999387/
      data (bm0ij (  9,  8,ibeta), ibeta = 1,10) /   &
        0.993250,  0.993837,  0.995007,  0.996208,  0.997223,   &
        0.998007,  0.998584,  0.998999,  0.999293,  0.999499/
      data (bm0ij (  9,  9,ibeta), ibeta = 1,10) /   &
        0.994413,  0.994903,  0.995878,  0.996876,  0.997716,   &
        0.998363,  0.998839,  0.999180,  0.999421,  0.999591/
      data (bm0ij (  9, 10,ibeta), ibeta = 1,10) /   &
        0.995376,  0.995785,  0.996597,  0.997425,  0.998121,   &
        0.998655,  0.999048,  0.999328,  0.999526,  0.999665/
      data (bm0ij ( 10,  1,ibeta), ibeta = 1,10) /   &
        0.989082,  0.989991,  0.991819,  0.993723,  0.995357,   &
        0.996637,  0.997592,  0.998286,  0.998781,  0.999132/
      data (bm0ij ( 10,  2,ibeta), ibeta = 1,10) /   &
        0.989613,  0.990480,  0.992224,  0.994039,  0.995594,   &
        0.996810,  0.997717,  0.998375,  0.998845,  0.999178/
      data (bm0ij ( 10,  3,ibeta), ibeta = 1,10) /   &
        0.990744,  0.991523,  0.993086,  0.994708,  0.996094,   &
        0.997176,  0.997981,  0.998564,  0.998980,  0.999274/
      data (bm0ij ( 10,  4,ibeta), ibeta = 1,10) /   &
        0.992041,  0.992716,  0.994070,  0.995470,  0.996662,   &
        0.997591,  0.998280,  0.998778,  0.999133,  0.999383/
      data (bm0ij ( 10,  5,ibeta), ibeta = 1,10) /   &
        0.993292,  0.993867,  0.995015,  0.996199,  0.997205,   &
        0.997985,  0.998564,  0.998981,  0.999277,  0.999487/
      data (bm0ij ( 10,  6,ibeta), ibeta = 1,10) /   &
        0.994411,  0.994894,  0.995857,  0.996847,  0.997685,   &
        0.998334,  0.998814,  0.999159,  0.999404,  0.999577/
      data (bm0ij ( 10,  7,ibeta), ibeta = 1,10) /   &
        0.995373,  0.995776,  0.996577,  0.997400,  0.998094,   &
        0.998630,  0.999026,  0.999310,  0.999512,  0.999654/
      data (bm0ij ( 10,  8,ibeta), ibeta = 1,10) /   &
        0.996181,  0.996516,  0.997181,  0.997861,  0.998435,   &
        0.998877,  0.999202,  0.999435,  0.999601,  0.999717/
      data (bm0ij ( 10,  9,ibeta), ibeta = 1,10) /   &
        0.996851,  0.997128,  0.997680,  0.998242,  0.998715,   &
        0.999079,  0.999346,  0.999538,  0.999673,  0.999769/
      data (bm0ij ( 10, 10,ibeta), ibeta = 1,10) /   &
        0.997402,  0.997632,  0.998089,  0.998554,  0.998945,   &
        0.999244,  0.999464,  0.999622,  0.999733,  0.999811/


! rpm....   3rd moment nuclei mode corr. fac. for bimodal fm coag rate

       data (bm3i( 1, 1,ibeta ), ibeta=1,10)/   &
       0.70708,0.71681,0.73821,0.76477,0.79350,0.82265,0.85090,0.87717,   &
       0.90069,0.92097/
       data (bm3i( 1, 2,ibeta ), ibeta=1,10)/   &
       0.72172,0.73022,0.74927,0.77324,0.79936,0.82601,0.85199,0.87637,   &
       0.89843,0.91774/
       data (bm3i( 1, 3,ibeta ), ibeta=1,10)/   &
       0.78291,0.78896,0.80286,0.82070,0.84022,0.85997,0.87901,0.89669,   &
       0.91258,0.92647/
       data (bm3i( 1, 4,ibeta ), ibeta=1,10)/   &
       0.87760,0.88147,0.89025,0.90127,0.91291,0.92420,0.93452,0.94355,   &
       0.95113,0.95726/
       data (bm3i( 1, 5,ibeta ), ibeta=1,10)/   &
       0.94988,0.95184,0.95612,0.96122,0.96628,0.97085,0.97467,0.97763,   &
       0.97971,0.98089/
       data (bm3i( 1, 6,ibeta ), ibeta=1,10)/   &
       0.98318,0.98393,0.98551,0.98728,0.98889,0.99014,0.99095,0.99124,   &
       0.99100,0.99020/
       data (bm3i( 1, 7,ibeta ), ibeta=1,10)/   &
       0.99480,0.99504,0.99551,0.99598,0.99629,0.99635,0.99611,0.99550,   &
       0.99450,0.99306/
       data (bm3i( 1, 8,ibeta ), ibeta=1,10)/   &
       0.99842,0.99848,0.99858,0.99861,0.99850,0.99819,0.99762,0.99674,   &
       0.99550,0.99388/
       data (bm3i( 1, 9,ibeta ), ibeta=1,10)/   &
       0.99951,0.99951,0.99949,0.99939,0.99915,0.99872,0.99805,0.99709,   &
       0.99579,0.99411/
       data (bm3i( 1,10,ibeta ), ibeta=1,10)/   &
       0.99984,0.99982,0.99976,0.99962,0.99934,0.99888,0.99818,0.99719,   &
       0.99587,0.99417/
       data (bm3i( 2, 1,ibeta ), ibeta=1,10)/   &
       0.72957,0.73993,0.76303,0.79178,0.82245,0.85270,0.88085,0.90578,   &
       0.92691,0.94415/
       data (bm3i( 2, 2,ibeta ), ibeta=1,10)/   &
       0.72319,0.73320,0.75547,0.78323,0.81307,0.84287,0.87107,0.89651,   &
       0.91852,0.93683/
       data (bm3i( 2, 3,ibeta ), ibeta=1,10)/   &
       0.74413,0.75205,0.76998,0.79269,0.81746,0.84258,0.86685,0.88938,   &
       0.90953,0.92695/
       data (bm3i( 2, 4,ibeta ), ibeta=1,10)/   &
       0.82588,0.83113,0.84309,0.85825,0.87456,0.89072,0.90594,0.91972,   &
       0.93178,0.94203/
       data (bm3i( 2, 5,ibeta ), ibeta=1,10)/   &
       0.91886,0.92179,0.92831,0.93624,0.94434,0.95192,0.95856,0.96409,   &
       0.96845,0.97164/
       data (bm3i( 2, 6,ibeta ), ibeta=1,10)/   &
       0.97129,0.97252,0.97515,0.97818,0.98108,0.98354,0.98542,0.98665,   &
       0.98721,0.98709/
       data (bm3i( 2, 7,ibeta ), ibeta=1,10)/   &
       0.99104,0.99145,0.99230,0.99320,0.99394,0.99439,0.99448,0.99416,   &
       0.99340,0.99217/
       data (bm3i( 2, 8,ibeta ), ibeta=1,10)/   &
       0.99730,0.99741,0.99763,0.99779,0.99782,0.99762,0.99715,0.99636,   &
       0.99519,0.99363/
       data (bm3i( 2, 9,ibeta ), ibeta=1,10)/   &
       0.99917,0.99919,0.99921,0.99915,0.99895,0.99856,0.99792,0.99698,   &
       0.99570,0.99404/
       data (bm3i( 2,10,ibeta ), ibeta=1,10)/   &
       0.99973,0.99973,0.99968,0.99955,0.99928,0.99883,0.99814,0.99716,   &
       0.99584,0.99415/
       data (bm3i( 3, 1,ibeta ), ibeta=1,10)/   &
       0.78358,0.79304,0.81445,0.84105,0.86873,0.89491,0.91805,0.93743,   &
       0.95300,0.96510/
       data (bm3i( 3, 2,ibeta ), ibeta=1,10)/   &
       0.76412,0.77404,0.79635,0.82404,0.85312,0.88101,0.90610,0.92751,   &
       0.94500,0.95879/
       data (bm3i( 3, 3,ibeta ), ibeta=1,10)/   &
       0.74239,0.75182,0.77301,0.79956,0.82809,0.85639,0.88291,0.90658,   &
       0.92683,0.94350/
       data (bm3i( 3, 4,ibeta ), ibeta=1,10)/   &
       0.78072,0.78758,0.80317,0.82293,0.84437,0.86589,0.88643,0.90526,   &
       0.92194,0.93625/
       data (bm3i( 3, 5,ibeta ), ibeta=1,10)/   &
       0.87627,0.88044,0.88981,0.90142,0.91357,0.92524,0.93585,0.94510,   &
       0.95285,0.95911/
       data (bm3i( 3, 6,ibeta ), ibeta=1,10)/   &
       0.95176,0.95371,0.95796,0.96297,0.96792,0.97233,0.97599,0.97880,   &
       0.98072,0.98178/
       data (bm3i( 3, 7,ibeta ), ibeta=1,10)/   &
       0.98453,0.98523,0.98670,0.98833,0.98980,0.99092,0.99160,0.99179,   &
       0.99145,0.99058/
       data (bm3i( 3, 8,ibeta ), ibeta=1,10)/   &
       0.99534,0.99555,0.99597,0.99637,0.99662,0.99663,0.99633,0.99569,   &
       0.99465,0.99318/
       data (bm3i( 3, 9,ibeta ), ibeta=1,10)/   &
       0.99859,0.99864,0.99872,0.99873,0.99860,0.99827,0.99768,0.99679,   &
       0.99555,0.99391/
       data (bm3i( 3,10,ibeta ), ibeta=1,10)/   &
       0.99956,0.99956,0.99953,0.99942,0.99918,0.99875,0.99807,0.99711,   &
       0.99580,0.99412/
       data (bm3i( 4, 1,ibeta ), ibeta=1,10)/   &
       0.84432,0.85223,0.86990,0.89131,0.91280,0.93223,0.94861,0.96172,   &
       0.97185,0.97945/
       data (bm3i( 4, 2,ibeta ), ibeta=1,10)/   &
       0.82299,0.83164,0.85101,0.87463,0.89857,0.92050,0.93923,0.95443,   &
       0.96629,0.97529/
       data (bm3i( 4, 3,ibeta ), ibeta=1,10)/   &
       0.77870,0.78840,0.81011,0.83690,0.86477,0.89124,0.91476,0.93460,   &
       0.95063,0.96316/
       data (bm3i( 4, 4,ibeta ), ibeta=1,10)/   &
       0.76386,0.77233,0.79147,0.81557,0.84149,0.86719,0.89126,0.91275,   &
       0.93116,0.94637/
       data (bm3i( 4, 5,ibeta ), ibeta=1,10)/   &
       0.82927,0.83488,0.84756,0.86346,0.88040,0.89704,0.91257,0.92649,   &
       0.93857,0.94874/
       data (bm3i( 4, 6,ibeta ), ibeta=1,10)/   &
       0.92184,0.92481,0.93136,0.93925,0.94724,0.95462,0.96104,0.96634,   &
       0.97048,0.97348/
       data (bm3i( 4, 7,ibeta ), ibeta=1,10)/   &
       0.97341,0.97457,0.97706,0.97991,0.98260,0.98485,0.98654,0.98760,   &
       0.98801,0.98777/
       data (bm3i( 4, 8,ibeta ), ibeta=1,10)/   &
       0.99192,0.99229,0.99305,0.99385,0.99449,0.99486,0.99487,0.99449,   &
       0.99367,0.99239/
       data (bm3i( 4, 9,ibeta ), ibeta=1,10)/   &
       0.99758,0.99768,0.99787,0.99800,0.99799,0.99777,0.99727,0.99645,   &
       0.99527,0.99369/
       data (bm3i( 4,10,ibeta ), ibeta=1,10)/   &
       0.99926,0.99928,0.99928,0.99921,0.99900,0.99860,0.99795,0.99701,   &
       0.99572,0.99405/
       data (bm3i( 5, 1,ibeta ), ibeta=1,10)/   &
       0.89577,0.90190,0.91522,0.93076,0.94575,0.95876,0.96932,0.97751,   &
       0.98367,0.98820/
       data (bm3i( 5, 2,ibeta ), ibeta=1,10)/   &
       0.87860,0.88547,0.90052,0.91828,0.93557,0.95075,0.96319,0.97292,   &
       0.98028,0.98572/
       data (bm3i( 5, 3,ibeta ), ibeta=1,10)/   &
       0.83381,0.84240,0.86141,0.88425,0.90707,0.92770,0.94510,0.95906,   &
       0.96986,0.97798/
       data (bm3i( 5, 4,ibeta ), ibeta=1,10)/   &
       0.78530,0.79463,0.81550,0.84127,0.86813,0.89367,0.91642,0.93566,   &
       0.95125,0.96347/
       data (bm3i( 5, 5,ibeta ), ibeta=1,10)/   &
       0.79614,0.80332,0.81957,0.84001,0.86190,0.88351,0.90368,0.92169,   &
       0.93718,0.95006/
       data (bm3i( 5, 6,ibeta ), ibeta=1,10)/   &
       0.88192,0.88617,0.89565,0.90728,0.91931,0.93076,0.94107,0.94997,   &
       0.95739,0.96333/
       data (bm3i( 5, 7,ibeta ), ibeta=1,10)/   &
       0.95509,0.95698,0.96105,0.96583,0.97048,0.97460,0.97796,0.98050,   &
       0.98218,0.98304/
       data (bm3i( 5, 8,ibeta ), ibeta=1,10)/   &
       0.98596,0.98660,0.98794,0.98943,0.99074,0.99172,0.99227,0.99235,   &
       0.99192,0.99096/
       data (bm3i( 5, 9,ibeta ), ibeta=1,10)/   &
       0.99581,0.99600,0.99637,0.99672,0.99691,0.99687,0.99653,0.99585,   &
       0.99478,0.99329/
       data (bm3i( 5,10,ibeta ), ibeta=1,10)/   &
       0.99873,0.99878,0.99884,0.99883,0.99869,0.99834,0.99774,0.99684,   &
       0.99558,0.99394/
       data (bm3i( 6, 1,ibeta ), ibeta=1,10)/   &
       0.93335,0.93777,0.94711,0.95764,0.96741,0.97562,0.98210,0.98701,   &
       0.99064,0.99327/
       data (bm3i( 6, 2,ibeta ), ibeta=1,10)/   &
       0.92142,0.92646,0.93723,0.94947,0.96096,0.97069,0.97842,0.98431,   &
       0.98868,0.99186/
       data (bm3i( 6, 3,ibeta ), ibeta=1,10)/   &
       0.88678,0.89351,0.90810,0.92508,0.94138,0.95549,0.96693,0.97578,   &
       0.98243,0.98731/
       data (bm3i( 6, 4,ibeta ), ibeta=1,10)/   &
       0.83249,0.84124,0.86051,0.88357,0.90655,0.92728,0.94477,0.95880,   &
       0.96964,0.97779/
       data (bm3i( 6, 5,ibeta ), ibeta=1,10)/   &
       0.79593,0.80444,0.82355,0.84725,0.87211,0.89593,0.91735,0.93566,   &
       0.95066,0.96255/
       data (bm3i( 6, 6,ibeta ), ibeta=1,10)/   &
       0.84124,0.84695,0.85980,0.87575,0.89256,0.90885,0.92383,0.93704,   &
       0.94830,0.95761/
       data (bm3i( 6, 7,ibeta ), ibeta=1,10)/   &
       0.92721,0.93011,0.93647,0.94406,0.95166,0.95862,0.96460,0.96949,   &
       0.97326,0.97595/
       data (bm3i( 6, 8,ibeta ), ibeta=1,10)/   &
       0.97573,0.97681,0.97913,0.98175,0.98421,0.98624,0.98772,0.98860,   &
       0.98885,0.98847/
       data (bm3i( 6, 9,ibeta ), ibeta=1,10)/   &
       0.99271,0.99304,0.99373,0.99444,0.99499,0.99528,0.99522,0.99477,   &
       0.99390,0.99258/
       data (bm3i( 6,10,ibeta ), ibeta=1,10)/   &
       0.99782,0.99791,0.99807,0.99817,0.99813,0.99788,0.99737,0.99653,   &
       0.99533,0.99374/
       data (bm3i( 7, 1,ibeta ), ibeta=1,10)/   &
       0.95858,0.96158,0.96780,0.97460,0.98073,0.98575,0.98963,0.99252,   &
       0.99463,0.99615/
       data (bm3i( 7, 2,ibeta ), ibeta=1,10)/   &
       0.95091,0.95438,0.96163,0.96962,0.97688,0.98286,0.98751,0.99099,   &
       0.99353,0.99536/
       data (bm3i( 7, 3,ibeta ), ibeta=1,10)/   &
       0.92751,0.93233,0.94255,0.95406,0.96473,0.97366,0.98070,0.98602,   &
       0.98994,0.99278/
       data (bm3i( 7, 4,ibeta ), ibeta=1,10)/   &
       0.88371,0.89075,0.90595,0.92351,0.94028,0.95474,0.96642,0.97544,   &
       0.98220,0.98715/
       data (bm3i( 7, 5,ibeta ), ibeta=1,10)/   &
       0.82880,0.83750,0.85671,0.87980,0.90297,0.92404,0.94195,0.95644,   &
       0.96772,0.97625/
       data (bm3i( 7, 6,ibeta ), ibeta=1,10)/   &
       0.81933,0.82655,0.84279,0.86295,0.88412,0.90449,0.92295,0.93890,   &
       0.95215,0.96281/
       data (bm3i( 7, 7,ibeta ), ibeta=1,10)/   &
       0.89099,0.89519,0.90448,0.91577,0.92732,0.93820,0.94789,0.95616,   &
       0.96297,0.96838/
       data (bm3i( 7, 8,ibeta ), ibeta=1,10)/   &
       0.95886,0.96064,0.96448,0.96894,0.97324,0.97701,0.98004,0.98228,   &
       0.98371,0.98435/
       data (bm3i( 7, 9,ibeta ), ibeta=1,10)/   &
       0.98727,0.98786,0.98908,0.99043,0.99160,0.99245,0.99288,0.99285,   &
       0.99234,0.99131/
       data (bm3i( 7,10,ibeta ), ibeta=1,10)/   &
       0.99621,0.99638,0.99671,0.99700,0.99715,0.99707,0.99670,0.99599,   &
       0.99489,0.99338/
       data (bm3i( 8, 1,ibeta ), ibeta=1,10)/   &
       0.97470,0.97666,0.98064,0.98491,0.98867,0.99169,0.99399,0.99569,   &
       0.99691,0.99779/
       data (bm3i( 8, 2,ibeta ), ibeta=1,10)/   &
       0.96996,0.97225,0.97693,0.98196,0.98643,0.99003,0.99279,0.99482,   &
       0.99630,0.99735/
       data (bm3i( 8, 3,ibeta ), ibeta=1,10)/   &
       0.95523,0.95848,0.96522,0.97260,0.97925,0.98468,0.98888,0.99200,   &
       0.99427,0.99590/
       data (bm3i( 8, 4,ibeta ), ibeta=1,10)/   &
       0.92524,0.93030,0.94098,0.95294,0.96397,0.97317,0.98038,0.98582,   &
       0.98981,0.99270/
       data (bm3i( 8, 5,ibeta ), ibeta=1,10)/   &
       0.87576,0.88323,0.89935,0.91799,0.93583,0.95126,0.96377,0.97345,   &
       0.98072,0.98606/
       data (bm3i( 8, 6,ibeta ), ibeta=1,10)/   &
       0.83078,0.83894,0.85705,0.87899,0.90126,0.92179,0.93950,0.95404,   &
       0.96551,0.97430/
       data (bm3i( 8, 7,ibeta ), ibeta=1,10)/   &
       0.85727,0.86294,0.87558,0.89111,0.90723,0.92260,0.93645,0.94841,   &
       0.95838,0.96643/
       data (bm3i( 8, 8,ibeta ), ibeta=1,10)/   &
       0.93337,0.93615,0.94220,0.94937,0.95647,0.96292,0.96840,0.97283,   &
       0.97619,0.97854/
       data (bm3i( 8, 9,ibeta ), ibeta=1,10)/   &
       0.97790,0.97891,0.98105,0.98346,0.98569,0.98751,0.98879,0.98950,   &
       0.98961,0.98912/
       data (bm3i( 8,10,ibeta ), ibeta=1,10)/   &
       0.99337,0.99367,0.99430,0.99493,0.99541,0.99562,0.99551,0.99501,   &
       0.99410,0.99274/
       data (bm3i( 9, 1,ibeta ), ibeta=1,10)/   &
       0.98470,0.98594,0.98844,0.99106,0.99334,0.99514,0.99650,0.99749,   &
       0.99821,0.99872/
       data (bm3i( 9, 2,ibeta ), ibeta=1,10)/   &
       0.98184,0.98330,0.98624,0.98934,0.99205,0.99420,0.99582,0.99701,   &
       0.99787,0.99848/
       data (bm3i( 9, 3,ibeta ), ibeta=1,10)/   &
       0.97288,0.97498,0.97927,0.98385,0.98789,0.99113,0.99360,0.99541,   &
       0.99673,0.99766/
       data (bm3i( 9, 4,ibeta ), ibeta=1,10)/   &
       0.95403,0.95741,0.96440,0.97202,0.97887,0.98444,0.98872,0.99190,   &
       0.99421,0.99586/
       data (bm3i( 9, 5,ibeta ), ibeta=1,10)/   &
       0.91845,0.92399,0.93567,0.94873,0.96076,0.97079,0.97865,0.98457,   &
       0.98892,0.99206/
       data (bm3i( 9, 6,ibeta ), ibeta=1,10)/   &
       0.86762,0.87533,0.89202,0.91148,0.93027,0.94669,0.96013,0.97062,   &
       0.97855,0.98441/
       data (bm3i( 9, 7,ibeta ), ibeta=1,10)/   &
       0.84550,0.85253,0.86816,0.88721,0.90671,0.92490,0.94083,0.95413,   &
       0.96481,0.97314/
       data (bm3i( 9, 8,ibeta ), ibeta=1,10)/   &
       0.90138,0.90544,0.91437,0.92513,0.93602,0.94615,0.95506,0.96258,   &
       0.96868,0.97347/
       data (bm3i( 9, 9,ibeta ), ibeta=1,10)/   &
       0.96248,0.96415,0.96773,0.97187,0.97583,0.97925,0.98198,0.98394,   &
       0.98514,0.98559/
       data (bm3i( 9,10,ibeta ), ibeta=1,10)/   &
       0.98837,0.98892,0.99005,0.99127,0.99232,0.99306,0.99339,0.99328,   &
       0.99269,0.99161/
       data (bm3i(10, 1,ibeta ), ibeta=1,10)/   &
       0.99080,0.99158,0.99311,0.99471,0.99607,0.99715,0.99795,0.99853,   &
       0.99895,0.99925/
       data (bm3i(10, 2,ibeta ), ibeta=1,10)/   &
       0.98910,0.99001,0.99182,0.99371,0.99533,0.99661,0.99757,0.99826,   &
       0.99876,0.99912/
       data (bm3i(10, 3,ibeta ), ibeta=1,10)/   &
       0.98374,0.98506,0.98772,0.99051,0.99294,0.99486,0.99630,0.99736,   &
       0.99812,0.99866/
       data (bm3i(10, 4,ibeta ), ibeta=1,10)/   &
       0.97238,0.97453,0.97892,0.98361,0.98773,0.99104,0.99354,0.99538,   &
       0.99671,0.99765/
       data (bm3i(10, 5,ibeta ), ibeta=1,10)/   &
       0.94961,0.95333,0.96103,0.96941,0.97693,0.98303,0.98772,0.99119,   &
       0.99371,0.99551/
       data (bm3i(10, 6,ibeta ), ibeta=1,10)/   &
       0.90943,0.91550,0.92834,0.94275,0.95608,0.96723,0.97600,0.98263,   &
       0.98751,0.99103/
       data (bm3i(10, 7,ibeta ), ibeta=1,10)/   &
       0.86454,0.87200,0.88829,0.90749,0.92630,0.94300,0.95687,0.96785,   &
       0.97626,0.98254/
       data (bm3i(10, 8,ibeta ), ibeta=1,10)/   &
       0.87498,0.88048,0.89264,0.90737,0.92240,0.93642,0.94877,0.95917,   &
       0.96762,0.97429/
       data (bm3i(10, 9,ibeta ), ibeta=1,10)/   &
       0.93946,0.94209,0.94781,0.95452,0.96111,0.96704,0.97203,0.97602,   &
       0.97900,0.98106/
       data (bm3i(10,10,ibeta ), ibeta=1,10)/   &
       0.97977,0.98071,0.98270,0.98492,0.98695,0.98858,0.98970,0.99027,   &
       0.99026,0.98968/

! fsb fm correction for intramodal m2 coagulation
       data bm2ii /   &
        0.707107,  0.720583,  0.745310,  0.748056,  0.696935,   &
        0.604164,  0.504622,  0.416559,  0.343394,  0.283641/

! *** total correction for intramodal m2 coagulation

      data bm2iitt /   &
        1.000000,  0.907452,  0.680931,  0.409815,  0.196425,   &
        0.078814,  0.028473,  0.009800,  0.003322,  0.001129/


! fsb fm correction for m2 i to j coagulation

      data (bm2ij (  1,  1,ibeta), ibeta = 1,10) /   &
        0.707107,  0.716828,  0.738240,  0.764827,  0.793610,   &
        0.822843,  0.851217,  0.877670,  0.901404,  0.921944/
      data (bm2ij (  1,  2,ibeta), ibeta = 1,10) /   &
        0.719180,  0.727975,  0.747638,  0.772334,  0.799234,   &
        0.826666,  0.853406,  0.878482,  0.901162,  0.920987/
      data (bm2ij (  1,  3,ibeta), ibeta = 1,10) /   &
        0.760947,  0.767874,  0.783692,  0.803890,  0.826015,   &
        0.848562,  0.870498,  0.891088,  0.909823,  0.926400/
      data (bm2ij (  1,  4,ibeta), ibeta = 1,10) /   &
        0.830926,  0.836034,  0.847708,  0.862528,  0.878521,   &
        0.894467,  0.909615,  0.923520,  0.935959,  0.946858/
      data (bm2ij (  1,  5,ibeta), ibeta = 1,10) /   &
        0.903643,  0.907035,  0.914641,  0.924017,  0.933795,   &
        0.943194,  0.951806,  0.959449,  0.966087,  0.971761/
      data (bm2ij (  1,  6,ibeta), ibeta = 1,10) /   &
        0.954216,  0.956094,  0.960211,  0.965123,  0.970068,   &
        0.974666,  0.978750,  0.982277,  0.985268,  0.987775/
      data (bm2ij (  1,  7,ibeta), ibeta = 1,10) /   &
        0.980546,  0.981433,  0.983343,  0.985568,  0.987751,   &
        0.989735,  0.991461,  0.992926,  0.994150,  0.995164/
      data (bm2ij (  1,  8,ibeta), ibeta = 1,10) /   &
        0.992142,  0.992524,  0.993338,  0.994272,  0.995174,   &
        0.995981,  0.996675,  0.997257,  0.997740,  0.998137/
      data (bm2ij (  1,  9,ibeta), ibeta = 1,10) /   &
        0.996868,  0.997026,  0.997361,  0.997742,  0.998106,   &
        0.998430,  0.998705,  0.998935,  0.999125,  0.999280/
      data (bm2ij (  1, 10,ibeta), ibeta = 1,10) /   &
        0.998737,  0.998802,  0.998939,  0.999094,  0.999241,   &
        0.999371,  0.999481,  0.999573,  0.999648,  0.999709/
      data (bm2ij (  2,  1,ibeta), ibeta = 1,10) /   &
        0.729600,  0.739948,  0.763059,  0.791817,  0.822510,   &
        0.852795,  0.881000,  0.905999,  0.927206,  0.944532/
      data (bm2ij (  2,  2,ibeta), ibeta = 1,10) /   &
        0.727025,  0.737116,  0.759615,  0.787657,  0.817740,   &
        0.847656,  0.875801,  0.901038,  0.922715,  0.940643/
      data (bm2ij (  2,  3,ibeta), ibeta = 1,10) /   &
        0.738035,  0.746779,  0.766484,  0.791340,  0.818324,   &
        0.845546,  0.871629,  0.895554,  0.916649,  0.934597/
      data (bm2ij (  2,  4,ibeta), ibeta = 1,10) /   &
        0.784185,  0.790883,  0.806132,  0.825501,  0.846545,   &
        0.867745,  0.888085,  0.906881,  0.923705,  0.938349/
      data (bm2ij (  2,  5,ibeta), ibeta = 1,10) /   &
        0.857879,  0.862591,  0.873238,  0.886539,  0.900645,   &
        0.914463,  0.927360,  0.939004,  0.949261,  0.958125/
      data (bm2ij (  2,  6,ibeta), ibeta = 1,10) /   &
        0.925441,  0.928304,  0.934645,  0.942324,  0.950181,   &
        0.957600,  0.964285,  0.970133,  0.975147,  0.979388/
      data (bm2ij (  2,  7,ibeta), ibeta = 1,10) /   &
        0.966728,  0.968176,  0.971323,  0.975027,  0.978705,   &
        0.982080,  0.985044,  0.987578,  0.989710,  0.991485/
      data (bm2ij (  2,  8,ibeta), ibeta = 1,10) /   &
        0.986335,  0.986980,  0.988362,  0.989958,  0.991511,   &
        0.992912,  0.994122,  0.995143,  0.995992,  0.996693/
      data (bm2ij (  2,  9,ibeta), ibeta = 1,10) /   &
        0.994547,  0.994817,  0.995391,  0.996046,  0.996677,   &
        0.997238,  0.997719,  0.998122,  0.998454,  0.998727/
      data (bm2ij (  2, 10,ibeta), ibeta = 1,10) /   &
        0.997817,  0.997928,  0.998163,  0.998429,  0.998683,   &
        0.998908,  0.999099,  0.999258,  0.999389,  0.999497/
      data (bm2ij (  3,  1,ibeta), ibeta = 1,10) /   &
        0.783612,  0.793055,  0.814468,  0.841073,  0.868769,   &
        0.894963,  0.918118,  0.937527,  0.953121,  0.965244/
      data (bm2ij (  3,  2,ibeta), ibeta = 1,10) /   &
        0.772083,  0.781870,  0.803911,  0.831238,  0.859802,   &
        0.887036,  0.911349,  0.931941,  0.948649,  0.961751/
      data (bm2ij (  3,  3,ibeta), ibeta = 1,10) /   &
        0.755766,  0.765509,  0.787380,  0.814630,  0.843526,   &
        0.871670,  0.897443,  0.919870,  0.938557,  0.953576/
      data (bm2ij (  3,  4,ibeta), ibeta = 1,10) /   &
        0.763816,  0.772145,  0.790997,  0.814784,  0.840434,   &
        0.865978,  0.890034,  0.911671,  0.930366,  0.945963/
      data (bm2ij (  3,  5,ibeta), ibeta = 1,10) /   &
        0.813597,  0.819809,  0.833889,  0.851618,  0.870640,   &
        0.889514,  0.907326,  0.923510,  0.937768,  0.950003/
      data (bm2ij (  3,  6,ibeta), ibeta = 1,10) /   &
        0.886317,  0.890437,  0.899643,  0.910955,  0.922730,   &
        0.934048,  0.944422,  0.953632,  0.961624,  0.968444/
      data (bm2ij (  3,  7,ibeta), ibeta = 1,10) /   &
        0.944565,  0.946855,  0.951872,  0.957854,  0.963873,   &
        0.969468,  0.974438,  0.978731,  0.982372,  0.985424/
      data (bm2ij (  3,  8,ibeta), ibeta = 1,10) /   &
        0.976358,  0.977435,  0.979759,  0.982467,  0.985125,   &
        0.987540,  0.989642,  0.991425,  0.992916,  0.994150/
      data (bm2ij (  3,  9,ibeta), ibeta = 1,10) /   &
        0.990471,  0.990932,  0.991917,  0.993048,  0.994142,   &
        0.995121,  0.995964,  0.996671,  0.997258,  0.997740/
      data (bm2ij (  3, 10,ibeta), ibeta = 1,10) /   &
        0.996199,  0.996389,  0.996794,  0.997254,  0.997694,   &
        0.998086,  0.998420,  0.998699,  0.998929,  0.999117/
      data (bm2ij (  4,  1,ibeta), ibeta = 1,10) /   &
        0.844355,  0.852251,  0.869914,  0.891330,  0.912823,   &
        0.932259,  0.948642,  0.961767,  0.971897,  0.979510/
      data (bm2ij (  4,  2,ibeta), ibeta = 1,10) /   &
        0.831550,  0.839954,  0.858754,  0.881583,  0.904592,   &
        0.925533,  0.943309,  0.957647,  0.968779,  0.977185/
      data (bm2ij (  4,  3,ibeta), ibeta = 1,10) /   &
        0.803981,  0.813288,  0.834060,  0.859400,  0.885285,   &
        0.909286,  0.930084,  0.947193,  0.960714,  0.971078/
      data (bm2ij (  4,  4,ibeta), ibeta = 1,10) /   &
        0.781787,  0.791080,  0.811931,  0.837749,  0.864768,   &
        0.890603,  0.913761,  0.933477,  0.949567,  0.962261/
      data (bm2ij (  4,  5,ibeta), ibeta = 1,10) /   &
        0.791591,  0.799355,  0.816916,  0.838961,  0.862492,   &
        0.885595,  0.907003,  0.925942,  0.942052,  0.955310/
      data (bm2ij (  4,  6,ibeta), ibeta = 1,10) /   &
        0.844933,  0.850499,  0.863022,  0.878593,  0.895038,   &
        0.911072,  0.925939,  0.939227,  0.950765,  0.960550/
      data (bm2ij (  4,  7,ibeta), ibeta = 1,10) /   &
        0.912591,  0.916022,  0.923607,  0.932777,  0.942151,   &
        0.951001,  0.958976,  0.965950,  0.971924,  0.976965/
      data (bm2ij (  4,  8,ibeta), ibeta = 1,10) /   &
        0.959859,  0.961617,  0.965433,  0.969924,  0.974382,   &
        0.978472,  0.982063,  0.985134,  0.987716,  0.989865/
      data (bm2ij (  4,  9,ibeta), ibeta = 1,10) /   &
        0.983377,  0.984162,  0.985844,  0.987788,  0.989681,   &
        0.991386,  0.992860,  0.994104,  0.995139,  0.995991/
      data (bm2ij (  4, 10,ibeta), ibeta = 1,10) /   &
        0.993343,  0.993672,  0.994370,  0.995169,  0.995937,   &
        0.996622,  0.997209,  0.997700,  0.998106,  0.998439/
      data (bm2ij (  5,  1,ibeta), ibeta = 1,10) /   &
        0.895806,  0.901918,  0.915233,  0.930783,  0.945768,   &
        0.958781,  0.969347,  0.977540,  0.983697,  0.988225/
      data (bm2ij (  5,  2,ibeta), ibeta = 1,10) /   &
        0.885634,  0.892221,  0.906629,  0.923540,  0.939918,   &
        0.954213,  0.965873,  0.974951,  0.981794,  0.986840/
      data (bm2ij (  5,  3,ibeta), ibeta = 1,10) /   &
        0.860120,  0.867858,  0.884865,  0.904996,  0.924724,   &
        0.942177,  0.956602,  0.967966,  0.976616,  0.983043/
      data (bm2ij (  5,  4,ibeta), ibeta = 1,10) /   &
        0.827462,  0.836317,  0.855885,  0.879377,  0.902897,   &
        0.924232,  0.942318,  0.956900,  0.968222,  0.976774/
      data (bm2ij (  5,  5,ibeta), ibeta = 1,10) /   &
        0.805527,  0.814279,  0.833853,  0.857892,  0.882726,   &
        0.906095,  0.926690,  0.943938,  0.957808,  0.968615/
      data (bm2ij (  5,  6,ibeta), ibeta = 1,10) /   &
        0.820143,  0.827223,  0.843166,  0.863002,  0.883905,   &
        0.904128,  0.922585,  0.938687,  0.952222,  0.963255/
      data (bm2ij (  5,  7,ibeta), ibeta = 1,10) /   &
        0.875399,  0.880208,  0.890929,  0.904065,  0.917699,   &
        0.930756,  0.942656,  0.953131,  0.962113,  0.969657/
      data (bm2ij (  5,  8,ibeta), ibeta = 1,10) /   &
        0.934782,  0.937520,  0.943515,  0.950656,  0.957840,   &
        0.964516,  0.970446,  0.975566,  0.979905,  0.983534/
      data (bm2ij (  5,  9,ibeta), ibeta = 1,10) /   &
        0.971369,  0.972679,  0.975505,  0.978797,  0.982029,   &
        0.984964,  0.987518,  0.989685,  0.991496,  0.992994/
      data (bm2ij (  5, 10,ibeta), ibeta = 1,10) /   &
        0.988329,  0.988893,  0.990099,  0.991485,  0.992825,   &
        0.994025,  0.995058,  0.995925,  0.996643,  0.997234/
      data (bm2ij (  6,  1,ibeta), ibeta = 1,10) /   &
        0.933384,  0.937784,  0.947130,  0.957655,  0.967430,   &
        0.975639,  0.982119,  0.987031,  0.990657,  0.993288/
      data (bm2ij (  6,  2,ibeta), ibeta = 1,10) /   &
        0.926445,  0.931227,  0.941426,  0.952975,  0.963754,   &
        0.972845,  0.980044,  0.985514,  0.989558,  0.992498/
      data (bm2ij (  6,  3,ibeta), ibeta = 1,10) /   &
        0.907835,  0.913621,  0.926064,  0.940308,  0.953745,   &
        0.965189,  0.974327,  0.981316,  0.986510,  0.990297/
      data (bm2ij (  6,  4,ibeta), ibeta = 1,10) /   &
        0.879088,  0.886306,  0.901945,  0.920079,  0.937460,   &
        0.952509,  0.964711,  0.974166,  0.981265,  0.986484/
      data (bm2ij (  6,  5,ibeta), ibeta = 1,10) /   &
        0.846500,  0.854862,  0.873189,  0.894891,  0.916264,   &
        0.935315,  0.951197,  0.963812,  0.973484,  0.980715/
      data (bm2ij (  6,  6,ibeta), ibeta = 1,10) /   &
        0.828137,  0.836250,  0.854310,  0.876287,  0.898710,   &
        0.919518,  0.937603,  0.952560,  0.964461,  0.973656/
      data (bm2ij (  6,  7,ibeta), ibeta = 1,10) /   &
        0.848595,  0.854886,  0.868957,  0.886262,  0.904241,   &
        0.921376,  0.936799,  0.950096,  0.961172,  0.970145/
      data (bm2ij (  6,  8,ibeta), ibeta = 1,10) /   &
        0.902919,  0.906922,  0.915760,  0.926427,  0.937312,   &
        0.947561,  0.956758,  0.964747,  0.971525,  0.977175/
      data (bm2ij (  6,  9,ibeta), ibeta = 1,10) /   &
        0.952320,  0.954434,  0.959021,  0.964418,  0.969774,   &
        0.974688,  0.979003,  0.982690,  0.985789,  0.988364/
      data (bm2ij (  6, 10,ibeta), ibeta = 1,10) /   &
        0.979689,  0.980650,  0.982712,  0.985093,  0.987413,   &
        0.989502,  0.991308,  0.992831,  0.994098,  0.995142/
      data (bm2ij (  7,  1,ibeta), ibeta = 1,10) /   &
        0.958611,  0.961598,  0.967817,  0.974620,  0.980752,   &
        0.985771,  0.989650,  0.992543,  0.994653,  0.996171/
      data (bm2ij (  7,  2,ibeta), ibeta = 1,10) /   &
        0.954225,  0.957488,  0.964305,  0.971795,  0.978576,   &
        0.984144,  0.988458,  0.991681,  0.994034,  0.995728/
      data (bm2ij (  7,  3,ibeta), ibeta = 1,10) /   &
        0.942147,  0.946158,  0.954599,  0.963967,  0.972529,   &
        0.979612,  0.985131,  0.989271,  0.992301,  0.994487/
      data (bm2ij (  7,  4,ibeta), ibeta = 1,10) /   &
        0.921821,  0.927048,  0.938140,  0.950598,  0.962118,   &
        0.971752,  0.979326,  0.985046,  0.989254,  0.992299/
      data (bm2ij (  7,  5,ibeta), ibeta = 1,10) /   &
        0.893419,  0.900158,  0.914598,  0.931070,  0.946584,   &
        0.959795,  0.970350,  0.978427,  0.984432,  0.988811/
      data (bm2ij (  7,  6,ibeta), ibeta = 1,10) /   &
        0.863302,  0.871111,  0.888103,  0.907990,  0.927305,   &
        0.944279,  0.958245,  0.969211,  0.977540,  0.983720/
      data (bm2ij (  7,  7,ibeta), ibeta = 1,10) /   &
        0.850182,  0.857560,  0.873890,  0.893568,  0.913408,   &
        0.931591,  0.947216,  0.960014,  0.970121,  0.977886/
      data (bm2ij (  7,  8,ibeta), ibeta = 1,10) /   &
        0.875837,  0.881265,  0.893310,  0.907936,  0.922910,   &
        0.936977,  0.949480,  0.960154,  0.968985,  0.976111/
      data (bm2ij (  7,  9,ibeta), ibeta = 1,10) /   &
        0.926228,  0.929445,  0.936486,  0.944868,  0.953293,   &
        0.961108,  0.968028,  0.973973,  0.978974,  0.983118/
      data (bm2ij (  7, 10,ibeta), ibeta = 1,10) /   &
        0.965533,  0.967125,  0.970558,  0.974557,  0.978484,   &
        0.982050,  0.985153,  0.987785,  0.989982,  0.991798/
      data (bm2ij (  8,  1,ibeta), ibeta = 1,10) /   &
        0.974731,  0.976674,  0.980660,  0.984926,  0.988689,   &
        0.991710,  0.994009,  0.995703,  0.996929,  0.997805/
      data (bm2ij (  8,  2,ibeta), ibeta = 1,10) /   &
        0.972062,  0.974192,  0.978571,  0.983273,  0.987432,   &
        0.990780,  0.993333,  0.995218,  0.996581,  0.997557/
      data (bm2ij (  8,  3,ibeta), ibeta = 1,10) /   &
        0.964662,  0.967300,  0.972755,  0.978659,  0.983921,   &
        0.988181,  0.991444,  0.993859,  0.995610,  0.996863/
      data (bm2ij (  8,  4,ibeta), ibeta = 1,10) /   &
        0.951782,  0.955284,  0.962581,  0.970559,  0.977737,   &
        0.983593,  0.988103,  0.991454,  0.993889,  0.995635/
      data (bm2ij (  8,  5,ibeta), ibeta = 1,10) /   &
        0.931947,  0.936723,  0.946751,  0.957843,  0.967942,   &
        0.976267,  0.982734,  0.987571,  0.991102,  0.993642/
      data (bm2ij (  8,  6,ibeta), ibeta = 1,10) /   &
        0.905410,  0.911665,  0.924950,  0.939908,  0.953798,   &
        0.965469,  0.974684,  0.981669,  0.986821,  0.990556/
      data (bm2ij (  8,  7,ibeta), ibeta = 1,10) /   &
        0.878941,  0.886132,  0.901679,  0.919688,  0.936970,   &
        0.951980,  0.964199,  0.973709,  0.980881,  0.986174/
      data (bm2ij (  8,  8,ibeta), ibeta = 1,10) /   &
        0.871653,  0.878218,  0.892652,  0.909871,  0.927034,   &
        0.942592,  0.955836,  0.966604,  0.975065,  0.981545/
      data (bm2ij (  8,  9,ibeta), ibeta = 1,10) /   &
        0.900693,  0.905239,  0.915242,  0.927232,  0.939335,   &
        0.950555,  0.960420,  0.968774,  0.975651,  0.981188/
      data (bm2ij (  8, 10,ibeta), ibeta = 1,10) /   &
        0.944922,  0.947435,  0.952894,  0.959317,  0.965689,   &
        0.971529,  0.976645,  0.981001,  0.984641,  0.987642/
      data (bm2ij (  9,  1,ibeta), ibeta = 1,10) /   &
        0.984736,  0.985963,  0.988453,  0.991078,  0.993357,   &
        0.995161,  0.996519,  0.997512,  0.998226,  0.998734/
      data (bm2ij (  9,  2,ibeta), ibeta = 1,10) /   &
        0.983141,  0.984488,  0.987227,  0.990119,  0.992636,   &
        0.994632,  0.996137,  0.997238,  0.998030,  0.998595/
      data (bm2ij (  9,  3,ibeta), ibeta = 1,10) /   &
        0.978726,  0.980401,  0.983819,  0.987450,  0.990626,   &
        0.993157,  0.995071,  0.996475,  0.997486,  0.998206/
      data (bm2ij (  9,  4,ibeta), ibeta = 1,10) /   &
        0.970986,  0.973224,  0.977818,  0.982737,  0.987072,   &
        0.990546,  0.993184,  0.995124,  0.996523,  0.997521/
      data (bm2ij (  9,  5,ibeta), ibeta = 1,10) /   &
        0.958579,  0.961700,  0.968149,  0.975116,  0.981307,   &
        0.986301,  0.990112,  0.992923,  0.994954,  0.996404/
      data (bm2ij (  9,  6,ibeta), ibeta = 1,10) /   &
        0.940111,  0.944479,  0.953572,  0.963506,  0.972436,   &
        0.979714,  0.985313,  0.989468,  0.992483,  0.994641/
      data (bm2ij (  9,  7,ibeta), ibeta = 1,10) /   &
        0.916127,  0.921878,  0.934003,  0.947506,  0.959899,   &
        0.970199,  0.978255,  0.984314,  0.988755,  0.991960/
      data (bm2ij (  9,  8,ibeta), ibeta = 1,10) /   &
        0.893848,  0.900364,  0.914368,  0.930438,  0.945700,   &
        0.958824,  0.969416,  0.977603,  0.983746,  0.988262/
      data (bm2ij (  9,  9,ibeta), ibeta = 1,10) /   &
        0.892161,  0.897863,  0.910315,  0.925021,  0.939523,   &
        0.952544,  0.963544,  0.972442,  0.979411,  0.984742/
      data (bm2ij (  9, 10,ibeta), ibeta = 1,10) /   &
        0.922260,  0.925966,  0.934047,  0.943616,  0.953152,   &
        0.961893,  0.969506,  0.975912,  0.981167,  0.985394/
      data (bm2ij ( 10,  1,ibeta), ibeta = 1,10) /   &
        0.990838,  0.991598,  0.993128,  0.994723,  0.996092,   &
        0.997167,  0.997969,  0.998552,  0.998969,  0.999265/
      data (bm2ij ( 10,  2,ibeta), ibeta = 1,10) /   &
        0.989892,  0.990727,  0.992411,  0.994167,  0.995678,   &
        0.996864,  0.997751,  0.998396,  0.998858,  0.999186/
      data (bm2ij ( 10,  3,ibeta), ibeta = 1,10) /   &
        0.987287,  0.988327,  0.990428,  0.992629,  0.994529,   &
        0.996026,  0.997148,  0.997965,  0.998551,  0.998967/
      data (bm2ij ( 10,  4,ibeta), ibeta = 1,10) /   &
        0.982740,  0.984130,  0.986952,  0.989926,  0.992508,   &
        0.994551,  0.996087,  0.997208,  0.998012,  0.998584/
      data (bm2ij ( 10,  5,ibeta), ibeta = 1,10) /   &
        0.975380,  0.977330,  0.981307,  0.985529,  0.989216,   &
        0.992147,  0.994358,  0.995975,  0.997136,  0.997961/
      data (bm2ij ( 10,  6,ibeta), ibeta = 1,10) /   &
        0.963911,  0.966714,  0.972465,  0.978614,  0.984022,   &
        0.988346,  0.991620,  0.994020,  0.995747,  0.996974/
      data (bm2ij ( 10,  7,ibeta), ibeta = 1,10) /   &
        0.947187,  0.951161,  0.959375,  0.968258,  0.976160,   &
        0.982540,  0.987409,  0.991000,  0.993592,  0.995441/
      data (bm2ij ( 10,  8,ibeta), ibeta = 1,10) /   &
        0.926045,  0.931270,  0.942218,  0.954297,  0.965273,   &
        0.974311,  0.981326,  0.986569,  0.990394,  0.993143/
      data (bm2ij ( 10,  9,ibeta), ibeta = 1,10) /   &
        0.908092,  0.913891,  0.926288,  0.940393,  0.953667,   &
        0.964987,  0.974061,  0.981038,  0.986253,  0.990078/
      data (bm2ij ( 10, 10,ibeta), ibeta = 1,10) /   &
        0.911143,  0.915972,  0.926455,  0.938721,  0.950701,   &
        0.961370,  0.970329,  0.977549,  0.983197,  0.987518/


! fsb total correction factor for m2 coagulation j from i

      data  (bm2ji( 1, 1,ibeta), ibeta = 1,10) /   &
        0.753466,  0.756888,  0.761008,  0.759432,  0.748675,   &
        0.726951,  0.693964,  0.650915,  0.600227,  0.545000/
      data  (bm2ji( 1, 2,ibeta), ibeta = 1,10) /   &
        0.824078,  0.828698,  0.835988,  0.838943,  0.833454,   &
        0.817148,  0.789149,  0.750088,  0.701887,  0.647308/
      data  (bm2ji( 1, 3,ibeta), ibeta = 1,10) /   &
        1.007389,  1.014362,  1.028151,  1.041011,  1.047939,   &
        1.045707,  1.032524,  1.007903,  0.972463,  0.927667/
      data  (bm2ji( 1, 4,ibeta), ibeta = 1,10) /   &
        1.246157,  1.255135,  1.274249,  1.295351,  1.313362,   &
        1.325187,  1.329136,  1.324491,  1.311164,  1.289459/
      data  (bm2ji( 1, 5,ibeta), ibeta = 1,10) /   &
        1.450823,  1.459551,  1.478182,  1.499143,  1.518224,   &
        1.533312,  1.543577,  1.548882,  1.549395,  1.545364/
      data  (bm2ji( 1, 6,ibeta), ibeta = 1,10) /   &
        1.575248,  1.581832,  1.595643,  1.610866,  1.624601,   &
        1.635690,  1.643913,  1.649470,  1.652688,  1.653878/
      data  (bm2ji( 1, 7,ibeta), ibeta = 1,10) /   &
        1.638426,  1.642626,  1.651293,  1.660641,  1.668926,   &
        1.675571,  1.680572,  1.684147,  1.686561,  1.688047/
      data  (bm2ji( 1, 8,ibeta), ibeta = 1,10) /   &
        1.669996,  1.672392,  1.677283,  1.682480,  1.687028,   &
        1.690651,  1.693384,  1.695372,  1.696776,  1.697734/
      data  (bm2ji( 1, 9,ibeta), ibeta = 1,10) /   &
        1.686148,  1.687419,  1.689993,  1.692704,  1.695057,   &
        1.696922,  1.698329,  1.699359,  1.700099,  1.700621/
      data  (bm2ji( 1,10,ibeta), ibeta = 1,10) /   &
        1.694364,  1.695010,  1.696313,  1.697676,  1.698853,   &
        1.699782,  1.700482,  1.700996,  1.701366,  1.701631/
      data  (bm2ji( 2, 1,ibeta), ibeta = 1,10) /   &
        0.783166,  0.779369,  0.768044,  0.747572,  0.716709,   &
        0.675422,  0.624981,  0.567811,  0.507057,  0.445975/
      data  (bm2ji( 2, 2,ibeta), ibeta = 1,10) /   &
        0.848390,  0.847100,  0.840874,  0.826065,  0.800296,   &
        0.762625,  0.713655,  0.655545,  0.591603,  0.525571/
      data  (bm2ji( 2, 3,ibeta), ibeta = 1,10) /   &
        1.039894,  1.043786,  1.049445,  1.049664,  1.039407,   &
        1.015322,  0.975983,  0.922180,  0.856713,  0.783634/
      data  (bm2ji( 2, 4,ibeta), ibeta = 1,10) /   &
        1.345995,  1.356064,  1.376947,  1.398304,  1.412685,   &
        1.414611,  1.400652,  1.369595,  1.322261,  1.260993/
      data  (bm2ji( 2, 5,ibeta), ibeta = 1,10) /   &
        1.675575,  1.689859,  1.720957,  1.756659,  1.788976,   &
        1.812679,  1.824773,  1.824024,  1.810412,  1.784630/
      data  (bm2ji( 2, 6,ibeta), ibeta = 1,10) /   &
        1.919835,  1.933483,  1.962973,  1.996810,  2.028377,   &
        2.054172,  2.072763,  2.083963,  2.088190,  2.086052/
      data  (bm2ji( 2, 7,ibeta), ibeta = 1,10) /   &
        2.064139,  2.074105,  2.095233,  2.118909,  2.140688,   &
        2.158661,  2.172373,  2.182087,  2.188330,  2.191650/
      data  (bm2ji( 2, 8,ibeta), ibeta = 1,10) /   &
        2.144871,  2.150990,  2.163748,  2.177731,  2.190364,   &
        2.200712,  2.208687,  2.214563,  2.218716,  2.221502/
      data  (bm2ji( 2, 9,ibeta), ibeta = 1,10) /   &
        2.189223,  2.192595,  2.199540,  2.207033,  2.213706,   &
        2.219125,  2.223297,  2.226403,  2.228660,  2.230265/
      data  (bm2ji( 2,10,ibeta), ibeta = 1,10) /   &
        2.212595,  2.214342,  2.217912,  2.221723,  2.225082,   &
        2.227791,  2.229869,  2.231417,  2.232551,  2.233372/
      data  (bm2ji( 3, 1,ibeta), ibeta = 1,10) /   &
        0.837870,  0.824476,  0.793119,  0.750739,  0.700950,   &
        0.646691,  0.590508,  0.534354,  0.479532,  0.426856/
      data  (bm2ji( 3, 2,ibeta), ibeta = 1,10) /   &
        0.896771,  0.885847,  0.859327,  0.821694,  0.775312,   &
        0.722402,  0.665196,  0.605731,  0.545742,  0.486687/
      data  (bm2ji( 3, 3,ibeta), ibeta = 1,10) /   &
        1.076089,  1.071727,  1.058845,  1.036171,  1.002539,   &
        0.957521,  0.901640,  0.836481,  0.764597,  0.689151/
      data  (bm2ji( 3, 4,ibeta), ibeta = 1,10) /   &
        1.409571,  1.415168,  1.425346,  1.432021,  1.428632,   &
        1.409696,  1.371485,  1.312958,  1.236092,  1.145293/
      data  (bm2ji( 3, 5,ibeta), ibeta = 1,10) /   &
        1.862757,  1.880031,  1.918394,  1.963456,  2.004070,   &
        2.030730,  2.036144,  2.016159,  1.970059,  1.900079/
      data  (bm2ji( 3, 6,ibeta), ibeta = 1,10) /   &
        2.289741,  2.313465,  2.366789,  2.431612,  2.495597,   &
        2.549838,  2.588523,  2.608665,  2.609488,  2.591662/
      data  (bm2ji( 3, 7,ibeta), ibeta = 1,10) /   &
        2.597157,  2.618731,  2.666255,  2.722597,  2.777531,   &
        2.825187,  2.862794,  2.889648,  2.906199,  2.913380/
      data  (bm2ji( 3, 8,ibeta), ibeta = 1,10) /   &
        2.797975,  2.813116,  2.845666,  2.882976,  2.918289,   &
        2.948461,  2.972524,  2.990687,  3.003664,  3.012284/
      data  (bm2ji( 3, 9,ibeta), ibeta = 1,10) /   &
        2.920832,  2.929843,  2.948848,  2.970057,  2.989632,   &
        3.006057,  3.019067,  3.028979,  3.036307,  3.041574/
      data  (bm2ji( 3,10,ibeta), ibeta = 1,10) /   &
        2.989627,  2.994491,  3.004620,  3.015720,  3.025789,   &
        3.034121,  3.040664,  3.045641,  3.049347,  3.052066/
      data  (bm2ji( 4, 1,ibeta), ibeta = 1,10) /   &
        0.893179,  0.870897,  0.820996,  0.759486,  0.695488,   &
        0.634582,  0.579818,  0.532143,  0.490927,  0.454618/
      data  (bm2ji( 4, 2,ibeta), ibeta = 1,10) /   &
        0.948355,  0.927427,  0.880215,  0.821146,  0.758524,   &
        0.697680,  0.641689,  0.591605,  0.546919,  0.506208/
      data  (bm2ji( 4, 3,ibeta), ibeta = 1,10) /   &
        1.109562,  1.093648,  1.056438,  1.007310,  0.951960,   &
        0.894453,  0.837364,  0.781742,  0.727415,  0.673614/
      data  (bm2ji( 4, 4,ibeta), ibeta = 1,10) /   &
        1.423321,  1.417557,  1.402442,  1.379079,  1.347687,   &
        1.308075,  1.259703,  1.201983,  1.134778,  1.058878/
      data  (bm2ji( 4, 5,ibeta), ibeta = 1,10) /   &
        1.933434,  1.944347,  1.968765,  1.997653,  2.023054,   &
        2.036554,  2.029949,  1.996982,  1.934982,  1.845473/
      data  (bm2ji( 4, 6,ibeta), ibeta = 1,10) /   &
        2.547772,  2.577105,  2.645918,  2.735407,  2.830691,   &
        2.917268,  2.981724,  3.013684,  3.007302,  2.961560/
      data  (bm2ji( 4, 7,ibeta), ibeta = 1,10) /   &
        3.101817,  3.139271,  3.225851,  3.336402,  3.453409,   &
        3.563116,  3.655406,  3.724014,  3.766113,  3.781394/
      data  (bm2ji( 4, 8,ibeta), ibeta = 1,10) /   &
        3.540920,  3.573780,  3.647439,  3.737365,  3.828468,   &
        3.911436,  3.981317,  4.036345,  4.076749,  4.103751/
      data  (bm2ji( 4, 9,ibeta), ibeta = 1,10) /   &
        3.856771,  3.879363,  3.928579,  3.986207,  4.042173,   &
        4.091411,  4.132041,  4.164052,  4.188343,  4.206118/
      data  (bm2ji( 4,10,ibeta), ibeta = 1,10) /   &
        4.053923,  4.067191,  4.095509,  4.127698,  4.158037,   &
        4.184055,  4.205135,  4.221592,  4.234115,  4.243463/
      data  (bm2ji( 5, 1,ibeta), ibeta = 1,10) /   &
        0.935846,  0.906814,  0.843358,  0.768710,  0.695885,   &
        0.631742,  0.579166,  0.538471,  0.508410,  0.486863/
      data  (bm2ji( 5, 2,ibeta), ibeta = 1,10) /   &
        0.988308,  0.959524,  0.896482,  0.821986,  0.748887,   &
        0.684168,  0.630908,  0.589516,  0.558676,  0.536056/
      data  (bm2ji( 5, 3,ibeta), ibeta = 1,10) /   &
        1.133795,  1.107139,  1.048168,  0.977258,  0.906341,   &
        0.842477,  0.789093,  0.746731,  0.713822,  0.687495/
      data  (bm2ji( 5, 4,ibeta), ibeta = 1,10) /   &
        1.405692,  1.385781,  1.340706,  1.284776,  1.227085,   &
        1.173532,  1.127008,  1.087509,  1.052712,  1.018960/
      data  (bm2ji( 5, 5,ibeta), ibeta = 1,10) /   &
        1.884992,  1.879859,  1.868463,  1.854995,  1.841946,   &
        1.829867,  1.816972,  1.799319,  1.771754,  1.729406/
      data  (bm2ji( 5, 6,ibeta), ibeta = 1,10) /   &
        2.592275,  2.612268,  2.661698,  2.731803,  2.815139,   &
        2.901659,  2.978389,  3.031259,  3.048045,  3.021122/
      data  (bm2ji( 5, 7,ibeta), ibeta = 1,10) /   &
        3.390321,  3.435519,  3.545615,  3.698419,  3.876958,   &
        4.062790,  4.236125,  4.378488,  4.475619,  4.519170/
      data  (bm2ji( 5, 8,ibeta), ibeta = 1,10) /   &
        4.161376,  4.216558,  4.346896,  4.519451,  4.711107,   &
        4.902416,  5.077701,  5.226048,  5.341423,  5.421764/
      data  (bm2ji( 5, 9,ibeta), ibeta = 1,10) /   &
        4.843961,  4.892035,  5.001492,  5.138515,  5.281684,   &
        5.416805,  5.535493,  5.634050,  5.712063,  5.770996/
      data  (bm2ji( 5,10,ibeta), ibeta = 1,10) /   &
        5.352093,  5.385119,  5.458056,  5.545311,  5.632162,   &
        5.710566,  5.777005,  5.830863,  5.873123,  5.905442/
      data  (bm2ji( 6, 1,ibeta), ibeta = 1,10) /   &
        0.964038,  0.930794,  0.859433,  0.777776,  0.700566,   &
        0.634671,  0.582396,  0.543656,  0.517284,  0.501694/
      data  (bm2ji( 6, 2,ibeta), ibeta = 1,10) /   &
        1.013416,  0.979685,  0.907197,  0.824135,  0.745552,   &
        0.678616,  0.625870,  0.587348,  0.561864,  0.547674/
      data  (bm2ji( 6, 3,ibeta), ibeta = 1,10) /   &
        1.145452,  1.111457,  1.038152,  0.953750,  0.873724,   &
        0.805955,  0.753621,  0.717052,  0.694920,  0.684910/
      data  (bm2ji( 6, 4,ibeta), ibeta = 1,10) /   &
        1.376547,  1.345004,  1.276415,  1.196704,  1.121091,   &
        1.058249,  1.012197,  0.983522,  0.970323,  0.968933/
      data  (bm2ji( 6, 5,ibeta), ibeta = 1,10) /   &
        1.778801,  1.755897,  1.706074,  1.649008,  1.597602,   &
        1.560087,  1.540365,  1.538205,  1.549738,  1.568333/
      data  (bm2ji( 6, 6,ibeta), ibeta = 1,10) /   &
        2.447603,  2.445172,  2.443762,  2.451842,  2.475877,   &
        2.519039,  2.580118,  2.653004,  2.727234,  2.789738/
      data  (bm2ji( 6, 7,ibeta), ibeta = 1,10) /   &
        3.368490,  3.399821,  3.481357,  3.606716,  3.772101,   &
        3.969416,  4.184167,  4.396163,  4.582502,  4.721838/
      data  (bm2ji( 6, 8,ibeta), ibeta = 1,10) /   &
        4.426458,  4.489861,  4.648250,  4.877510,  5.160698,   &
        5.477495,  5.803123,  6.111250,  6.378153,  6.586050/
      data  (bm2ji( 6, 9,ibeta), ibeta = 1,10) /   &
        5.568061,  5.644988,  5.829837,  6.081532,  6.371214,   &
        6.672902,  6.963737,  7.226172,  7.449199,  7.627886/
      data  (bm2ji( 6,10,ibeta), ibeta = 1,10) /   &
        6.639152,  6.707020,  6.863974,  7.065285,  7.281744,   &
        7.492437,  7.683587,  7.847917,  7.983296,  8.090977/
      data  (bm2ji( 7, 1,ibeta), ibeta = 1,10) /   &
        0.980853,  0.945724,  0.871244,  0.787311,  0.708818,   &
        0.641987,  0.588462,  0.547823,  0.518976,  0.500801/
      data  (bm2ji( 7, 2,ibeta), ibeta = 1,10) /   &
        1.026738,  0.990726,  0.914306,  0.828140,  0.747637,   &
        0.679351,  0.625127,  0.584662,  0.556910,  0.540749/
      data  (bm2ji( 7, 3,ibeta), ibeta = 1,10) /   &
        1.146496,  1.108808,  1.028695,  0.938291,  0.854101,   &
        0.783521,  0.728985,  0.690539,  0.667272,  0.657977/
      data  (bm2ji( 7, 4,ibeta), ibeta = 1,10) /   &
        1.344846,  1.306434,  1.224543,  1.132031,  1.046571,   &
        0.976882,  0.926488,  0.896067,  0.884808,  0.891027/
      data  (bm2ji( 7, 5,ibeta), ibeta = 1,10) /   &
        1.670227,  1.634583,  1.558421,  1.472939,  1.396496,   &
        1.339523,  1.307151,  1.300882,  1.319622,  1.360166/
      data  (bm2ji( 7, 6,ibeta), ibeta = 1,10) /   &
        2.224548,  2.199698,  2.148284,  2.095736,  2.059319,   &
        2.050496,  2.075654,  2.136382,  2.229641,  2.347958/
      data  (bm2ji( 7, 7,ibeta), ibeta = 1,10) /   &
        3.104483,  3.105947,  3.118398,  3.155809,  3.230427,   &
        3.350585,  3.519071,  3.731744,  3.976847,  4.235616/
      data  (bm2ji( 7, 8,ibeta), ibeta = 1,10) /   &
        4.288426,  4.331456,  4.447024,  4.633023,  4.891991,   &
        5.221458,  5.610060,  6.036467,  6.471113,  6.880462/
      data  (bm2ji( 7, 9,ibeta), ibeta = 1,10) /   &
        5.753934,  5.837061,  6.048530,  6.363800,  6.768061,   &
        7.241280,  7.755346,  8.276666,  8.771411,  9.210826/
      data  (bm2ji( 7,10,ibeta), ibeta = 1,10) /   &
        7.466219,  7.568810,  7.819032,  8.168340,  8.582973,   &
        9.030174,  9.478159,  9.899834, 10.275940, 10.595910/
      data  (bm2ji( 8, 1,ibeta), ibeta = 1,10) /   &
        0.990036,  0.954782,  0.880531,  0.797334,  0.719410,   &
        0.652220,  0.596923,  0.552910,  0.519101,  0.494529/
      data  (bm2ji( 8, 2,ibeta), ibeta = 1,10) /   &
        1.032428,  0.996125,  0.919613,  0.833853,  0.753611,   &
        0.684644,  0.628260,  0.583924,  0.550611,  0.527407/
      data  (bm2ji( 8, 3,ibeta), ibeta = 1,10) /   &
        1.141145,  1.102521,  1.021017,  0.929667,  0.844515,   &
        0.772075,  0.714086,  0.670280,  0.639824,  0.621970/
      data  (bm2ji( 8, 4,ibeta), ibeta = 1,10) /   &
        1.314164,  1.273087,  1.186318,  1.089208,  0.999476,   &
        0.924856,  0.867948,  0.829085,  0.807854,  0.803759/
      data  (bm2ji( 8, 5,ibeta), ibeta = 1,10) /   &
        1.580611,  1.538518,  1.449529,  1.350459,  1.260910,   &
        1.190526,  1.143502,  1.121328,  1.124274,  1.151974/
      data  (bm2ji( 8, 6,ibeta), ibeta = 1,10) /   &
        2.016773,  1.977721,  1.895727,  1.806974,  1.732891,   &
        1.685937,  1.673026,  1.697656,  1.761039,  1.862391/
      data  (bm2ji( 8, 7,ibeta), ibeta = 1,10) /   &
        2.750093,  2.723940,  2.672854,  2.628264,  2.612250,   &
        2.640406,  2.723211,  2.866599,  3.071893,  3.335217/
      data  (bm2ji( 8, 8,ibeta), ibeta = 1,10) /   &
        3.881905,  3.887143,  3.913667,  3.981912,  4.111099,   &
        4.316575,  4.608146,  4.988157,  5.449592,  5.974848/
      data  (bm2ji( 8, 9,ibeta), ibeta = 1,10) /   &
        5.438870,  5.492742,  5.640910,  5.886999,  6.241641,   &
        6.710609,  7.289480,  7.960725,  8.693495,  9.446644/
      data  (bm2ji( 8,10,ibeta), ibeta = 1,10) /   &
        7.521152,  7.624621,  7.892039,  8.300444,  8.839787,   &
        9.493227, 10.231770, 11.015642, 11.799990, 12.542260/
      data  (bm2ji( 9, 1,ibeta), ibeta = 1,10) /   &
        0.994285,  0.960012,  0.887939,  0.807040,  0.730578,   &
        0.663410,  0.606466,  0.559137,  0.520426,  0.489429/
      data  (bm2ji( 9, 2,ibeta), ibeta = 1,10) /   &
        1.033505,  0.998153,  0.923772,  0.840261,  0.761383,   &
        0.692242,  0.633873,  0.585709,  0.546777,  0.516215/
      data  (bm2ji( 9, 3,ibeta), ibeta = 1,10) /   &
        1.132774,  1.094907,  1.015161,  0.925627,  0.841293,   &
        0.767888,  0.706741,  0.657439,  0.619135,  0.591119/
      data  (bm2ji( 9, 4,ibeta), ibeta = 1,10) /   &
        1.286308,  1.245273,  1.158809,  1.061889,  0.971208,   &
        0.893476,  0.830599,  0.782561,  0.748870,  0.729198/
      data  (bm2ji( 9, 5,ibeta), ibeta = 1,10) /   &
        1.511105,  1.467141,  1.374520,  1.271162,  1.175871,   &
        1.096887,  1.037243,  0.997820,  0.978924,  0.980962/
      data  (bm2ji( 9, 6,ibeta), ibeta = 1,10) /   &
        1.857468,  1.812177,  1.717002,  1.612197,  1.519171,   &
        1.448660,  1.405871,  1.393541,  1.413549,  1.467532/
      data  (bm2ji( 9, 7,ibeta), ibeta = 1,10) /   &
        2.430619,  2.388452,  2.301326,  2.210241,  2.139724,   &
        2.104571,  2.114085,  2.174696,  2.291294,  2.467500/
      data  (bm2ji( 9, 8,ibeta), ibeta = 1,10) /   &
        3.385332,  3.357690,  3.306611,  3.269804,  3.274462,   &
        3.340862,  3.484609,  3.717740,  4.048748,  4.481588/
      data  (bm2ji( 9, 9,ibeta), ibeta = 1,10) /   &
        4.850497,  4.858280,  4.896008,  4.991467,  5.171511,   &
        5.459421,  5.873700,  6.426128,  7.119061,  7.942603/
      data  (bm2ji( 9,10,ibeta), ibeta = 1,10) /   &
        6.957098,  7.020164,  7.197272,  7.499331,  7.946554,   &
        8.555048,  9.330503, 10.263610, 11.327454, 12.478332/
      data  (bm2ji(10, 1,ibeta), ibeta = 1,10) /   &
        0.994567,  0.961842,  0.892854,  0.814874,  0.740198,   &
        0.673303,  0.615105,  0.565139,  0.522558,  0.486556/
      data  (bm2ji(10, 2,ibeta), ibeta = 1,10) /   &
        1.031058,  0.997292,  0.926082,  0.845571,  0.768501,   &
        0.699549,  0.639710,  0.588538,  0.545197,  0.508894/
      data  (bm2ji(10, 3,ibeta), ibeta = 1,10) /   &
        1.122535,  1.086287,  1.009790,  0.923292,  0.840626,   &
        0.766982,  0.703562,  0.650004,  0.605525,  0.569411/
      data  (bm2ji(10, 4,ibeta), ibeta = 1,10) /   &
        1.261142,  1.221555,  1.137979,  1.043576,  0.953745,   &
        0.874456,  0.807292,  0.752109,  0.708326,  0.675477/
      data  (bm2ji(10, 5,ibeta), ibeta = 1,10) /   &
        1.456711,  1.413432,  1.322096,  1.219264,  1.122319,   &
        1.038381,  0.969743,  0.916811,  0.879544,  0.858099/
      data  (bm2ji(10, 6,ibeta), ibeta = 1,10) /   &
        1.741792,  1.695157,  1.596897,  1.487124,  1.385734,   &
        1.301670,  1.238638,  1.198284,  1.181809,  1.190689/
      data  (bm2ji(10, 7,ibeta), ibeta = 1,10) /   &
        2.190197,  2.141721,  2.040226,  1.929245,  1.832051,   &
        1.760702,  1.721723,  1.719436,  1.757705,  1.840677/
      data  (bm2ji(10, 8,ibeta), ibeta = 1,10) /   &
        2.940764,  2.895085,  2.801873,  2.707112,  2.638603,   &
        2.613764,  2.644686,  2.741255,  2.912790,  3.168519/
      data  (bm2ji(10, 9,ibeta), ibeta = 1,10) /   &
        4.186191,  4.155844,  4.101953,  4.069102,  4.089886,   &
        4.189530,  4.389145,  4.707528,  5.161567,  5.765283/
      data  (bm2ji(10,10,ibeta), ibeta = 1,10) /   &
        6.119526,  6.127611,  6.171174,  6.286528,  6.508738,   &
        6.869521,  7.396912,  8.113749,  9.034683, 10.162190/

! *** end of data statements.


! *** start calculations:

      constii = abs( half * ( two ) ** two3rds - one )
      sqrttwo = sqrt(two)
      dlgsqt2 = one / log( sqrttwo )

         esat01   = exp( 0.125d0 * xxlsgat * xxlsgat )
         esac01   = exp( 0.125d0 * xxlsgac * xxlsgac )

         esat04  = esat01 ** 4
         esac04  = esac01 ** 4

         esat05  = esat04 * esat01
         esac05  = esac04 * esac01

         esat08  = esat04 * esat04
         esac08  = esac04 * esac04

         esat09  = esat08 * esat01
         esac09  = esac08 * esac01

         esat16  = esat08 * esat08
         esac16  = esac08 * esac08

         esat20  = esat16 * esat04
         esac20  = esac16 * esac04

         esat24  = esat20 * esat04
         esac24  = esac20 * esac04

         esat25  = esat20 * esat05
         esac25  = esac20 * esac05

         esat36  = esat20 * esat16
         esac36  = esac20 * esac16

         esat49  = esat24 * esat25

         esat64  = esat20 * esat20 * esat24
         esac64  = esac20 * esac20 * esac24

         esat100 = esat64 * esat36

         dgat2   = dgatk * dgatk
         dgat3   = dgatk * dgatk * dgatk
         dgac2   = dgacc * dgacc
         dgac3   = dgacc * dgacc * dgacc

         sqdgat  = sqrt( dgatk )
         sqdgac  = sqrt( dgacc )
         sqdgat5 = dgat2 * sqdgat
         sqdgac5 = dgac2 * sqdgac
         sqdgat7 = dgat3 * sqdgat

         xm2at = dgat2 * esat16
         xm3at = dgat3 * esat36

         xm2ac = dgac2 * esac16
         xm3ac = dgac3 * esac36

! *** for the free molecular regime:  page h.3 of whitby et al. (1991)

         r       = sqdgac / sqdgat
         r2      = r * r
         r3      = r2 * r
         rx4     = r2 * r2
         r5      = r3 * r2
         r6      = r3 * r3
         rx8      = rx4 * rx4
         ri1     = one / r
         ri2     = one / r2
         ri3     = one / r3
         ri4     = ri2 * ri2
         kngat   = two * lamda / dgatk
         kngac   = two * lamda / dgacc


! *** calculate ratio of geometric mean diameters
         rat = dgacc / dgatk
! *** trap subscripts for bm0 and bm0i, between 1 and 10
!     see page h.5 of whitby et al. (1991)

      n2n = max( 1, min( 10,   &
            nint( 4.0 * ( sgatk - 0.75d0 ) ) ) )

      n2a = max( 1, min( 10,   &
            nint( 4.0 * ( sgacc - 0.75d0 ) ) ) )

      n1  = max( 1, min( 10,   &
             1 + nint( dlgsqt2 * log( rat ) ) ) )

! *** intermodal coagulation


! *** set up for zeroeth moment

! *** near-continuum form:  equation h.10a of whitby et al. (1991)

         coagnc0 = knc * (   &
          two + a * ( kngat * ( esat04 + r2 * esat16 * esac04 )   &
                    + kngac * ( esac04 + ri2 * esac16 * esat04 ) )   &
                    + ( r2 + ri2 ) * esat04 * esac04  )


! *** free-molecular form:  equation h.7a of whitby et al. (1991)

         coagfm0 = kfmatac * sqdgat * bm0ij(n1,n2n,n2a) * (   &
                   esat01 + r * esac01 + two * r2 * esat01 * esac04   &
                 + rx4 * esat09 * esac16 + ri3 * esat16 * esac09   &
                 + two * ri1 * esat04 + esac01  )


! *** loss to accumulation mode

! *** harmonic mean

      coagatac0 = coagnc0 * coagfm0 / ( coagnc0 + coagfm0 )

      qn12 = coagatac0


! *** set up for second moment
!      the second moment equations are new and begin with equations a1
!     through a4 of binkowski and shankar (1995). after some algebraic
!     rearrangement and application of the extended mean value theorem
!     of integral calculus, equations are obtained that can be solved
!     analytically with correction factors as has been done by
!     whitby et al. (1991)

! *** the term ( dp1 + dp2 ) ** (2/3) in equations a3 and a4 of
!     binkowski and shankar (1995) is approximated by
!     (dgat ** 3 + dgac **3 ) ** 2/3

! *** near-continuum form

      i1nc = knc * dgat2 * (   &
             two * esat16   &
           + r2 * esat04 * esac04   &
           + ri2 * esat36 * esac04   &
           + a * kngat * (   &
                 esat04   &
           +     ri2 * esat16 * esac04   &
           +     ri4 * esat36 * esac16   &
           +     r2 * esac04 )  )




! *** free-molecular form

       i1fm =  kfmatac * sqdgat5 * bm2ij(n1,n2n,n2a) * (   &
               esat25   &
            +  two * r2 * esat09 * esac04   &
            +  rx4 * esat01 * esac16   &
            +  ri3 * esat64 * esac09   &
            +  two * ri1 * esat36 * esac01   &
            +  r * esat16 * esac01  )



! *** loss to accumulation mode

! *** harmonic mean

      i1 = ( i1fm * i1nc ) / ( i1fm + i1nc )

      coagatac2 = i1

      qs12 = coagatac2


! *** gain by accumulation mode

      coagacat2 = ( ( one + r6 ) ** two3rds - rx4 ) * i1

      qs21 = coagacat2 * bm2ji(n1,n2n,n2a)

! *** set up for third moment

! *** near-continuum form: equation h.10b of whitby et al. (1991)

      coagnc3 = knc * dgat3 * (   &
                two * esat36   &
              + a * kngat * ( esat16 + r2 * esat04 * esac04 )   &
              + a * kngac * ( esat36 * esac04 + ri2 * esat64 * esac16 )   &
              + r2 * esat16 * esac04 + ri2 * esat64 * esac04 )


! *** free_molecular form: equation h.7b of whitby et al. (1991)

      coagfm3 = kfmatac * sqdgat7 * bm3i( n1, n2n, n2a ) * (   &
               esat49   &
              +  r * esat36  * esac01   &
              + two * r2 * esat25  * esac04   &
              + rx4 * esat09  * esac16   &
              + ri3 * esat100 * esac09   &
              + two * ri1 * esat64  * esac01 )

! *** gain by accumulation mode = loss from aitken mode

! *** harmonic mean

      coagatac3 = coagnc3 * coagfm3 / ( coagnc3 + coagfm3 )

      qv12 = coagatac3

! *** intramodal coagulation

! *** zeroeth moment

! *** aitken mode

! *** near-continuum form: equation h.12a of whitby et al. (1991)

      coagnc_at = knc * (one + esat08 + a * kngat * (esat20 + esat04))

! *** free-molecular form: equation h.11a of whitby et al. (1991)

      coagfm_at = kfmat * sqdgat * bm0(n2n) *   &
                 ( esat01 + esat25 + two * esat05 )


! *** harmonic mean

      coagatat0 = coagfm_at * coagnc_at / ( coagfm_at + coagnc_at )

      qn11 = coagatat0


! *** accumulation mode

! *** near-continuum form: equation h.12a of whitby et al. (1991)

      coagnc_ac = knc * (one + esac08 + a * kngac * (esac20 + esac04))

! *** free-molecular form: equation h.11a of whitby et al. (1991)

      coagfm_ac = kfmac * sqdgac * bm0(n2a) *   &
                   ( esac01 + esac25 + two * esac05 )

! *** harmonic mean

      coagacac0 = coagfm_ac * coagnc_ac / ( coagfm_ac + coagnc_ac )

      qn22 = coagacac0


! *** set up for second moment
!      the second moment equations are new and begin with 3.11a on page
!     45 of whitby et al. (1991). after some algebraic rearrangement and
!     application of the extended mean value theorem of integral calculus
!     equations are obtained that can be solved analytically with
!     correction factors as has been done by whitby et al. (1991)

! *** aitken mode

! *** near-continuum

      i1nc_at = knc * dgat2 * (   &
             two * esat16   &
           + esat04 * esat04   &
           + esat36 * esat04   &
           + a * kngat * (   &
                two * esat04   &
           +     esat16 * esat04   &
           +     esat36 * esat16 )  )

! *** free- molecular form

       i1fm_at =  kfmat * sqdgat5 * bm2ii(n2n) * (   &
               esat25   &
            +  two * esat09 * esat04   &
            +  esat01 * esat16   &
            +  esat64 * esat09   &
            +  two * esat36 * esat01   &
            +  esat16 * esat01  )

      i1_at = ( i1nc_at * i1fm_at ) / ( i1nc_at + i1fm_at  )

      coagatat2 = constii * i1_at

      qs11 = coagatat2 * bm2iitt(n2n)

! *** accumulation mode

! *** near-continuum

      i1nc_ac = knc * dgac2 * (   &
             two * esac16   &
           + esac04 * esac04   &
           + esac36 * esac04   &
           + a * kngac * (   &
                two * esac04   &
           +     esac16 * esac04   &
           +     esac36 * esac16 )  )

! *** free- molecular form

       i1fm_ac =  kfmac * sqdgac5 * bm2ii(n2a) * (   &
               esac25   &
            +  two * esac09 * esac04   &
            +  esac01 * esac16   &
            +  esac64 * esac09   &
            +  two * esac36 * esac01   &
            +  esac16 * esac01  )

      i1_ac = ( i1nc_ac * i1fm_ac ) / ( i1nc_ac + i1fm_ac  )

      coagacac2 = constii * i1_ac

      qs22 = coagacac2 * bm2iitt(n2a)


      return

      end  subroutine getcoags

!----------------------------------------------------------------------
!----------------------------------------------------------------------

   end module modal_aero_coag



